Madrid House Price
Sergio Fernández, Miguel Ocón, Enrique Roa
26 de abril de 2023
1 Introducción y definición de objetivos
El análisis de datos en el mercado inmobiliario es algo muy común y que lleva muchísimos años desarrollándose. Con el objetivo de predecir que aspectos influyen principalmente en el precio de las casas en Madrid, hemos seleccionado un dataset (kaggle - Madrid House Price) con viviendas en venta de la capital española.
library(readr)
library(dplyr)
library(tidyr)
library(kableExtra)
library(ggplot2)
library(GGally)
library(gridExtra)
library(cowplot)
library(ggcorrplot)
library(gmodels)
library(caret)
library(ggfortify)
library(scales)
library(class)
library(distances)
library(visreg)
library(rpart)
library(rpart.plot)
library(rattle)
library(randomForest)
library(e1071)
library(kernlab)
library(pROC)
library(cluster)
library(tidyverse)
library(stringr)
library(Metrics)
library(factoextra)
library(NbClust)paste(R.Version()$version.string)## [1] "R version 4.1.2 (2021-11-01)"2 Análisis exploratorio inicial
2.1 Lectura y preparación de los datos
En esta primera etapa cargamos el dataset descargado y hacemos una observación superficial de las variables recogidas y sus características. Se representan también las 10 filas iniciales.
mhp <- read_csv("house_price_madrid_14_08_2022.csv")head(mhp, 10) %>%
kbl() %>%
kable_material(c("striped", "hover")) %>%
scroll_box(width = "100%", height = "350px")| price | house_type | house_type_2 | rooms | m2 | elevator | garage | neighborhood | district |
|---|---|---|---|---|---|---|---|---|
| 495000 | planta 1 | exterior | 3 | 118 | TRUE | TRUE | Chopera | Arganzuela |
| 485000 | planta 2 | exterior | 2 | 82 | TRUE | TRUE | Palos de Moguer | Arganzuela |
| 315000 | planta 2 | exterior | 2 | 72 | FALSE | FALSE | Legazpi | Arganzuela |
| 585000 | planta 4 | exterior | 2 | 174 | TRUE | TRUE | Palos de Moguer | Arganzuela |
| 255000 | bajo | exterior | 3 | 75 | FALSE | FALSE | Acacias | Arganzuela |
| 299000 | planta 1 | exterior | 1 | 69 | TRUE | FALSE | Chopera | Arganzuela |
| 265000 | planta 3 | exterior | 2 | 54 | TRUE | FALSE | Delicias | Arganzuela |
| 290000 | planta 3 | interior | 4 | 69 | TRUE | FALSE | Chopera | Arganzuela |
| 660220 | planta 2 | exterior | 3 | 129 | TRUE | TRUE | Imperial | Arganzuela |
| 525000 | planta 4 | exterior | 4 | 111 | TRUE | FALSE | Chopera | Arganzuela |
summary(mhp)## price house_type house_type_2 rooms
## Min. : 725 Length:15975 Length:15975 Min. : 1.000
## 1st Qu.: 195000 Class :character Class :character 1st Qu.: 2.000
## Median : 359973 Mode :character Mode :character Median : 3.000
## Mean : 624233 Mean : 2.848
## 3rd Qu.: 749000 3rd Qu.: 3.000
## Max. :13950000 Max. :41.000
## m2 elevator garage neighborhood
## Min. : 1.0 Mode :logical Mode :logical Length:15975
## 1st Qu.: 66.0 FALSE:4597 FALSE:11528 Class :character
## Median : 93.0 TRUE :11378 TRUE :4447 Mode :character
## Mean :124.8
## 3rd Qu.:142.0
## Max. :989.0
## district
## Length:15975
## Class :character
## Mode :character
##
##
## str(mhp) ## spc_tbl_ [15,975 × 9] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ price : num [1:15975] 495000 485000 315000 585000 255000 ...
## $ house_type : chr [1:15975] "planta 1" "planta 2" "planta 2" "planta 4" ...
## $ house_type_2: chr [1:15975] "exterior" "exterior" "exterior" "exterior" ...
## $ rooms : num [1:15975] 3 2 2 2 3 1 2 4 3 4 ...
## $ m2 : num [1:15975] 118 82 72 174 75 69 54 69 129 111 ...
## $ elevator : logi [1:15975] TRUE TRUE FALSE TRUE FALSE TRUE ...
## $ garage : logi [1:15975] TRUE TRUE FALSE TRUE FALSE FALSE ...
## $ neighborhood: chr [1:15975] "Chopera" "Palos de Moguer" "Legazpi" "Palos de Moguer" ...
## $ district : chr [1:15975] "Arganzuela" "Arganzuela" "Arganzuela" "Arganzuela" ...
## - attr(*, "spec")=
## .. cols(
## .. price = col_double(),
## .. house_type = col_character(),
## .. house_type_2 = col_character(),
## .. rooms = col_double(),
## .. m2 = col_double(),
## .. elevator = col_logical(),
## .. garage = col_logical(),
## .. neighborhood = col_character(),
## .. district = col_character()
## .. )
## - attr(*, "problems")=<externalptr>El dataset contiene 15.975 observaciones (correspondientes a una vivienda cada una) y 9 variables (de las cuales son 6 cualitativas y 3 cuantitativas).
A continuación, la descripción de cada una de las variables:
price: precio
house_type: tipo de vivienda (casa, chalet, piso…)
house_type_2: si es exterior o interior
rooms: número de habitaciones
m2: metros cuadrados
elevator: si tiene ascensor
garage: si incluye garaje
neighborhood: barrio de Madrid
district: distrito de Madrid
Tras esta primera observación se harán algunas modificaciones en el dataset para transformar a variables de tipo categórico las que correspondan, así como transformar a binarias (1/0) aquellas con sólo dos valores posibles.
# Preparación de los datos. Convertimos a factores las variables que son categóricas y a 1-0 las binarias.
mhp <- mhp %>% rename(exterior = house_type_2)
mhp$exterior = ifelse(mhp$exterior == "exterior", 1, 0)
mhp$exterior = factor(mhp$exterior, levels = c(1,0))
mhp$elevator = ifelse(mhp$elevator == "TRUE", 1, 0)
mhp$elevator = factor(mhp$elevator, levels = c(1,0))
mhp$garage = ifelse(mhp$garage == "TRUE", 1, 0)
mhp$garage = factor(mhp$garage, levels = c(1,0))
mhp$house_type <- as.factor(mhp$house_type)
mhp$exterior <- as.factor(mhp$exterior)
mhp$elevator <- as.factor(mhp$elevator)
mhp$garage <- as.factor(mhp$garage)
mhp$neighborhood <- as.factor(mhp$neighborhood)
mhp$district <- as.factor(mhp$district)2.2 Tratamiento de datos faltantes
Convertimos en NA aquellas observaciones que claramente contienen algún error:
La casa de 41 habitaciones (probablemente fueran 4 y esté mal imputado)
Las que tienen menos de 3 m2 (probablemente mal imputados al usar el punto para separar los miles)
La casa con un precio de 725 (claramente equivocado).
mhp$rooms[mhp$rooms > 40] <- NA
mhp$m2[mhp$m2 < 3] <- NA
mhp$price[mhp$price < 1000] <- NAPor otro lado, la variable house_type_2 tiene 469 filas sin datos. Junto a los NA añadidos en el paso anterior son el 3.2% del total, por lo que optaremos por eliminar estas filas del dataset.
sum(is.na(mhp))/nrow(mhp)## [1] 0.03205008mhp <- na.omit(mhp)2.3 División del dataset
Para finalizar la preparación de los datos, dividiremos nuestro dataset en tres partes, una de entrenamiento (train), con el que trabajaremos durante todo el análisis, otra de prueba (test) que nos servirá para comprobar la eficacia de los diferentes modelos predictivos que se apliquen, y finalmente uno de validación (validation), que se mantendrá sin usar hasta el último día del proyecto y con el cual validaremos el modelo final.
set.seed(108)
numero_total = nrow(mhp)
# Porcentajes de train, test y validation
w_train = .5
w_test = .25
w_validation = 1 - (w_train + w_test)
# Todos los índices
indices = seq(1:numero_total)
# Muestreo
indices_train = sample(1:numero_total, numero_total * w_train)
indices_test = sample(indices[-indices_train], numero_total * w_test)
indices_validation = indices[-c(indices_train,indices_test)]
# Agrupamos
mhp_train = mhp[indices_train,]
mhp_test = mhp[indices_test,]
mhp_validation = mhp[indices_validation,]3 Análisis univariante
3.1 Análisis de variables cualitativas
Una vez preparado el dataset comenzamos con el análisis de las variables cualitativas, que son las siguientes:
Tipo de vivienda
merge(setNames(as.data.frame(table(mhp_train$house_type)), c("house_type", "count")),
setNames(as.data.frame(round(prop.table(table(mhp_train$house_type))*100, 2)), c("house_type", "prop (%)"))
) %>%
arrange(desc(count)) %>%
kbl() %>%
kable_material(c("striped", "hover")) %>%
scroll_box(width = "100%", height = "350px")| house_type | count | prop (%) |
|---|---|---|
| planta 1 | 1581 | 20.45 |
| planta 2 | 1332 | 17.23 |
| planta 3 | 1180 | 15.26 |
| bajo | 996 | 12.88 |
| planta 4 | 801 | 10.36 |
| planta 5 | 545 | 7.05 |
| planta 6 | 317 | 4.10 |
| planta 7 | 196 | 2.53 |
| chalet | 180 | 2.33 |
| entreplanta | 131 | 1.69 |
| planta 8 | 119 | 1.54 |
| planta 9 | 78 | 1.01 |
| casa | 65 | 0.84 |
| semi-sotano | 48 | 0.62 |
| planta 10 | 44 | 0.57 |
| planta 12 | 23 | 0.30 |
| planta 13 | 22 | 0.28 |
| planta 11 | 21 | 0.27 |
| sotano | 19 | 0.25 |
| planta 14 | 10 | 0.13 |
| planta -1 | 8 | 0.10 |
| planta 15 | 6 | 0.08 |
| planta 16 | 4 | 0.05 |
| planta 18 | 2 | 0.03 |
| planta 19 | 2 | 0.03 |
| planta 20 | 2 | 0.03 |
ggplot(mhp_train, aes(house_type)) +
geom_bar(fill = "#0BB363") +
coord_flip() +
labs(x = "Tipo de vivienda", y = "Número de viviendas", title = "Viviendas por tipo") +
theme(plot.title = element_text(hjust = 0.5))
Como era de esperar, las viviendas más habituales son pisos en plantas no demasiado altas (bajos, primeros, segundos…). El número de casas y chalets también es pequeño en proporción.
Orientación exterior
merge(setNames(as.data.frame(table(mhp_train$exterior)), c("exterior", "count")),
setNames(as.data.frame(round(prop.table(table(mhp_train$exterior))*100, 2)), c("exterior", "prop (%)"))
) %>%
kbl() %>%
kable_material(c("striped", "hover"))| exterior | count | prop (%) |
|---|---|---|
| 0 | 848 | 10.97 |
| 1 | 6884 | 89.03 |
ggplot(mhp_train, aes(exterior)) +
geom_bar(fill = "#0BB363") +
labs(x = "Exterior", y = "Número de vivienda de viviendas", title = "Viviendas exteriores") +
theme(plot.title = element_text(hjust = 0.5))
Son una inmensa mayoría los pisos que dan a exterior.
Ascensor
merge(setNames(as.data.frame(table(mhp_train$elevator)), c("elevator", "count")),
setNames(as.data.frame(round(prop.table(table(mhp_train$elevator))*100, 2)), c("elevator", "prop (%)"))
) %>%
kbl() %>%
kable_material(c("striped", "hover"))| elevator | count | prop (%) |
|---|---|---|
| 0 | 2169 | 28.05 |
| 1 | 5563 | 71.95 |
ggplot(mhp_train, aes(elevator)) +
geom_bar(fill = "#0BB363") +
labs(x = "Ascensor", y = "Número de viviendas", title = "Viviendas con ascensor") +
theme(plot.title = element_text(hjust = 0.5))
También es más frecuente que tengan ascensor, aunque casi el 30% no lo tienen.
Garaje
merge(setNames(as.data.frame(table(mhp_train$garage)), c("garage", "count")),
setNames(as.data.frame(round(prop.table(table(mhp_train$garage))*100, 2)), c("garage", "prop (%)"))
) %>%
kbl() %>%
kable_material(c("striped", "hover"))| garage | count | prop (%) |
|---|---|---|
| 0 | 5570 | 72.04 |
| 1 | 2162 | 27.96 |
ggplot(mhp_train, aes(garage)) +
geom_bar(fill = "#0BB363") +
labs(x = "Garaje", y = "Número de viviendas", title = "Viviendas con garaje") +
theme(plot.title = element_text(hjust = 0.5))
Parece normal que la mayoría de pisos no tuvieran garaje incluido.
Barrio
Esta variable no será muy útil para el análisis, ya que en muchos casos son descripciones de la vivienda hechas por el usuario, por lo que habría que hacer un tratamiento previo para poder agruparlas. Algo innecesario ya que ese mismo efecto se consigue con el siguiente campo: distritos.
merge(setNames(as.data.frame(table(mhp_train$neighborhood)), c("neighborhood", "count")),
setNames(as.data.frame(round(prop.table(table(mhp_train$neighborhood))*100, 2)), c("neighborhood", "prop (%)"))
) %>%
kbl() %>%
kable_material(c("striped", "hover")) %>%
scroll_box(width = "100%", height = "350px")| neighborhood | count | prop (%) |
|---|---|---|
| 12 de Octubre-Orcasur | 19 | 0.25 |
| Abrantes | 54 | 0.70 |
| Acacias | 50 | 0.65 |
| Adelfas | 30 | 0.39 |
| Aeropuerto | 2 | 0.03 |
| Águilas | 33 | 0.43 |
| Alameda de Osuna | 32 | 0.41 |
| Almagro | 140 | 1.81 |
| Almendrales | 42 | 0.54 |
| Aluche | 56 | 0.72 |
| Ambroz | 30 | 0.39 |
| Amposta | 18 | 0.23 |
| Apóstol Santiago | 19 | 0.25 |
| Arapiles | 67 | 0.87 |
| Aravaca | 47 | 0.61 |
| Arcos | 26 | 0.34 |
| Argüelles | 142 | 1.84 |
| Arroyo del Fresno | 11 | 0.14 |
| Atalaya | 6 | 0.08 |
| Ático en Almagro | 8 | 0.10 |
| Ático en Arapiles | 1 | 0.01 |
| Ático en Aravaca | 0 | 0.00 |
| Ático en Arcos | 1 | 0.01 |
| Ático en Argüelles | 0 | 0.00 |
| Ático en Bernabéu-Hispanoamérica | 1 | 0.01 |
| Ático en Campo de las Naciones-Corralejos | 1 | 0.01 |
| Ático en Castellana | 10 | 0.13 |
| Ático en Castilla | 4 | 0.05 |
| Ático en Chueca-Justicia | 1 | 0.01 |
| Ático en Ciudad Universitaria | 1 | 0.01 |
| Ático en Colina | 1 | 0.01 |
| Ático en Concepción | 1 | 0.01 |
| Ático en Conde Orgaz-Piovera | 1 | 0.01 |
| Ático en Costillares | 1 | 0.01 |
| Ático en Cuatro Caminos | 0 | 0.00 |
| Ático en Delicias | 1 | 0.01 |
| Ático en El Cañaveral | 1 | 0.01 |
| Ático en El Viso | 8 | 0.10 |
| Ático en Fuente del Berro | 0 | 0.00 |
| Ático en Gaztambide | 2 | 0.03 |
| Ático en Goya | 3 | 0.04 |
| Ático en Guindalera | 2 | 0.03 |
| Ático en Horcajo | 1 | 0.01 |
| Ático en Ibiza | 2 | 0.03 |
| Ático en Jerónimos | 1 | 0.01 |
| Ático en Las Tablas | 1 | 0.01 |
| Ático en Lista | 1 | 0.01 |
| Ático en Malasaña-Universidad | 1 | 0.01 |
| Ático en Niño Jesús | 1 | 0.01 |
| Ático en Nueva España | 4 | 0.05 |
| Ático en Nuevos Ministerios-Ríos Rosas | 3 | 0.04 |
| Ático en Pacífico | 1 | 0.01 |
| Ático en Palacio | 0 | 0.00 |
| Ático en Palomas | 1 | 0.01 |
| Ático en Peñagrande | 1 | 0.01 |
| Ático en Pinar del Rey | 0 | 0.00 |
| Ático en Pueblo Nuevo | 1 | 0.01 |
| Ático en Puerta Bonita | 0 | 0.00 |
| Ático en Puerta del Ángel | 1 | 0.01 |
| Ático en Recoletos | 3 | 0.04 |
| Ático en Rejas | 0 | 0.00 |
| Ático en Rosas | 1 | 0.01 |
| Ático en Salvador | 1 | 0.01 |
| Ático en San Juan Bautista | 1 | 0.01 |
| Ático en San Pascual | 1 | 0.01 |
| Ático en Simancas | 0 | 0.00 |
| Ático en Sol | 3 | 0.04 |
| Ático en Trafalgar | 1 | 0.01 |
| Ático en Valdeacederas | 1 | 0.01 |
| Ático en Valdebebas - Valdefuentes | 1 | 0.01 |
| Ático en Valdebernardo - Valderrivas | 0 | 0.00 |
| Ático en Valdemarín | 2 | 0.03 |
| Ático en Vallehermoso | 2 | 0.03 |
| Bellas Vistas | 89 | 1.15 |
| Bernabéu-Hispanoamérica | 91 | 1.18 |
| Berruguete | 82 | 1.06 |
| Buena Vista | 40 | 0.52 |
| Butarque | 14 | 0.18 |
| Campamento | 16 | 0.21 |
| Campo de las Naciones-Corralejos | 16 | 0.21 |
| Canillas | 43 | 0.56 |
| Canillejas | 78 | 1.01 |
| Casa de Campo | 18 | 0.23 |
| Casa o chalet independiente en Aravaca | 7 | 0.09 |
| Casa o chalet independiente en Bernabéu-Hispanoamérica | 2 | 0.03 |
| Casa o chalet independiente en Canillas | 2 | 0.03 |
| Casa o chalet independiente en Canillejas | 1 | 0.01 |
| Casa o chalet independiente en Casa de Campo | 0 | 0.00 |
| Casa o chalet independiente en Casco Histórico de Barajas | 0 | 0.00 |
| Casa o chalet independiente en Castilla | 0 | 0.00 |
| Casa o chalet independiente en Ciudad Universitaria | 5 | 0.06 |
| Casa o chalet independiente en Conde Orgaz-Piovera | 7 | 0.09 |
| Casa o chalet independiente en El Plantío | 1 | 0.01 |
| Casa o chalet independiente en El Viso | 4 | 0.05 |
| Casa o chalet independiente en Fuentelarreina | 0 | 0.00 |
| Casa o chalet independiente en Guindalera | 0 | 0.00 |
| Casa o chalet independiente en Mirasierra | 4 | 0.05 |
| Casa o chalet independiente en Niño Jesús | 1 | 0.01 |
| Casa o chalet independiente en Nueva España | 6 | 0.08 |
| Casa o chalet independiente en Opañel | 0 | 0.00 |
| Casa o chalet independiente en Orcasitas | 0 | 0.00 |
| Casa o chalet independiente en Peñagrande | 4 | 0.05 |
| Casa o chalet independiente en Rejas | 1 | 0.01 |
| Casa o chalet independiente en Tres Olivos - Valverde | 2 | 0.03 |
| Casa o chalet independiente en Valdemarín | 7 | 0.09 |
| Casa o chalet independiente en Valdezarza | 0 | 0.00 |
| Casa o chalet independiente en Villaverde Alto | 0 | 0.00 |
| Casa o chalet independiente en Vista Alegre | 0 | 0.00 |
| Casco Histórico de Barajas | 22 | 0.28 |
| Casco Histórico de Vallecas | 87 | 1.13 |
| Casco Histórico de Vicálvaro | 35 | 0.45 |
| Castellana | 172 | 2.22 |
| Castilla | 62 | 0.80 |
| Chalet adosado en Alameda de Osuna | 1 | 0.01 |
| Chalet adosado en Aravaca | 0 | 0.00 |
| Chalet adosado en Arroyo del Fresno | 1 | 0.01 |
| Chalet adosado en Bernabéu-Hispanoamérica | 1 | 0.01 |
| Chalet adosado en Buena Vista | 1 | 0.01 |
| Chalet adosado en Campo de las Naciones-Corralejos | 1 | 0.01 |
| Chalet adosado en Canillas | 1 | 0.01 |
| Chalet adosado en Canillejas | 0 | 0.00 |
| Chalet adosado en Ciudad Jardín | 1 | 0.01 |
| Chalet adosado en Ciudad Universitaria | 0 | 0.00 |
| Chalet adosado en Concepción | 1 | 0.01 |
| Chalet adosado en Conde Orgaz-Piovera | 5 | 0.06 |
| Chalet adosado en El Plantío | 2 | 0.03 |
| Chalet adosado en El Viso | 4 | 0.05 |
| Chalet adosado en Fuente del Berro | 1 | 0.01 |
| Chalet adosado en Guindalera | 0 | 0.00 |
| Chalet adosado en Los Cármenes | 0 | 0.00 |
| Chalet adosado en Mirasierra | 1 | 0.01 |
| Chalet adosado en Moscardó | 1 | 0.01 |
| Chalet adosado en Nueva España | 3 | 0.04 |
| Chalet adosado en Numancia | 0 | 0.00 |
| Chalet adosado en Orcasitas | 2 | 0.03 |
| Chalet adosado en Palomas | 1 | 0.01 |
| Chalet adosado en Peñagrande | 1 | 0.01 |
| Chalet adosado en Pinar del Rey | 1 | 0.01 |
| Chalet adosado en Puerta Bonita | 0 | 0.00 |
| Chalet adosado en Rejas | 1 | 0.01 |
| Chalet adosado en Salvador | 0 | 0.00 |
| Chalet adosado en Tres Olivos - Valverde | 0 | 0.00 |
| Chalet adosado en Valdeacederas | 0 | 0.00 |
| Chalet adosado en Valdebebas - Valdefuentes | 1 | 0.01 |
| Chalet adosado en Valdezarza | 1 | 0.01 |
| Chalet en Águilas | 0 | 0.00 |
| Chalet en Aravaca | 1 | 0.01 |
| Chalet en Ciudad Universitaria | 1 | 0.01 |
| Chalet en El Plantío | 1 | 0.01 |
| Chalet en El Viso | 0 | 0.00 |
| Chalet en Entrevías | 1 | 0.01 |
| Chalet en Fuentelarreina | 0 | 0.00 |
| Chalet en Mirasierra | 1 | 0.01 |
| Chalet en Rejas | 0 | 0.00 |
| Chalet en Valdemarín | 0 | 0.00 |
| Chalet pareado en Abrantes | 1 | 0.01 |
| Chalet pareado en Alameda de Osuna | 0 | 0.00 |
| Chalet pareado en Apóstol Santiago | 1 | 0.01 |
| Chalet pareado en Aravaca | 3 | 0.04 |
| Chalet pareado en Arroyo del Fresno | 1 | 0.01 |
| Chalet pareado en Bernabéu-Hispanoamérica | 0 | 0.00 |
| Chalet pareado en Campo de las Naciones-Corralejos | 0 | 0.00 |
| Chalet pareado en Canillas | 0 | 0.00 |
| Chalet pareado en Canillejas | 0 | 0.00 |
| Chalet pareado en Ciudad Jardín | 0 | 0.00 |
| Chalet pareado en Ciudad Universitaria | 1 | 0.01 |
| Chalet pareado en Colina | 1 | 0.01 |
| Chalet pareado en Conde Orgaz-Piovera | 4 | 0.05 |
| Chalet pareado en El Plantío | 2 | 0.03 |
| Chalet pareado en El Viso | 3 | 0.04 |
| Chalet pareado en Fuentelarreina | 0 | 0.00 |
| Chalet pareado en Guindalera | 0 | 0.00 |
| Chalet pareado en Mirasierra | 3 | 0.04 |
| Chalet pareado en Nueva España | 0 | 0.00 |
| Chalet pareado en Palomas | 4 | 0.05 |
| Chalet pareado en Peñagrande | 1 | 0.01 |
| Chalet pareado en Rejas | 1 | 0.01 |
| Chalet pareado en Tres Olivos - Valverde | 0 | 0.00 |
| Chalet pareado en Valdebebas - Valdefuentes | 1 | 0.01 |
| Chopera | 34 | 0.44 |
| Chueca-Justicia | 67 | 0.87 |
| Ciudad Jardín | 47 | 0.61 |
| Ciudad Universitaria | 40 | 0.52 |
| Colina | 25 | 0.32 |
| Comillas | 52 | 0.67 |
| Concepción | 49 | 0.63 |
| Conde Orgaz-Piovera | 35 | 0.45 |
| Costillares | 53 | 0.69 |
| Cuatro Caminos | 98 | 1.27 |
| Cuatro Vientos | 1 | 0.01 |
| Cuzco-Castillejos | 98 | 1.27 |
| Delicias | 58 | 0.75 |
| Dúplex en Almagro | 6 | 0.08 |
| Dúplex en Apóstol Santiago | 0 | 0.00 |
| Dúplex en Arapiles | 0 | 0.00 |
| Dúplex en Aravaca | 3 | 0.04 |
| Dúplex en Arcos | 0 | 0.00 |
| Dúplex en Argüelles | 0 | 0.00 |
| Dúplex en Arroyo del Fresno | 1 | 0.01 |
| Dúplex en Bernabéu-Hispanoamérica | 0 | 0.00 |
| Dúplex en Berruguete | 1 | 0.01 |
| Dúplex en Butarque | 0 | 0.00 |
| Dúplex en Canillas | 1 | 0.01 |
| Dúplex en Canillejas | 1 | 0.01 |
| Dúplex en Casco Histórico de Barajas | 0 | 0.00 |
| Dúplex en Casco Histórico de Vallecas | 0 | 0.00 |
| Dúplex en Castellana | 1 | 0.01 |
| Dúplex en Castilla | 1 | 0.01 |
| Dúplex en Chueca-Justicia | 1 | 0.01 |
| Dúplex en Colina | 1 | 0.01 |
| Dúplex en Concepción | 2 | 0.03 |
| Dúplex en Conde Orgaz-Piovera | 1 | 0.01 |
| Dúplex en Costillares | 0 | 0.00 |
| Dúplex en Cuzco-Castillejos | 1 | 0.01 |
| Dúplex en El Viso | 4 | 0.05 |
| Dúplex en Entrevías | 1 | 0.01 |
| Dúplex en Fontarrón | 1 | 0.01 |
| Dúplex en Fuente del Berro | 0 | 0.00 |
| Dúplex en Goya | 1 | 0.01 |
| Dúplex en Imperial | 1 | 0.01 |
| Dúplex en Jerónimos | 4 | 0.05 |
| Dúplex en Las Tablas | 2 | 0.03 |
| Dúplex en Legazpi | 2 | 0.03 |
| Dúplex en Lista | 0 | 0.00 |
| Dúplex en Malasaña-Universidad | 0 | 0.00 |
| Dúplex en Mirasierra | 1 | 0.01 |
| Dúplex en Moscardó | 1 | 0.01 |
| Dúplex en Nueva España | 8 | 0.10 |
| Dúplex en Nuevos Ministerios-Ríos Rosas | 1 | 0.01 |
| Dúplex en Numancia | 2 | 0.03 |
| Dúplex en Opañel | 1 | 0.01 |
| Dúplex en Pacífico | 3 | 0.04 |
| Dúplex en Palacio | 1 | 0.01 |
| Dúplex en Pinar del Rey | 0 | 0.00 |
| Dúplex en Prosperidad | 1 | 0.01 |
| Dúplex en Puerta Bonita | 0 | 0.00 |
| Dúplex en Puerta del Ángel | 0 | 0.00 |
| Dúplex en Quintana | 0 | 0.00 |
| Dúplex en Recoletos | 3 | 0.04 |
| Dúplex en Rejas | 0 | 0.00 |
| Dúplex en San Isidro | 1 | 0.01 |
| Dúplex en San Juan Bautista | 6 | 0.08 |
| Dúplex en Sanchinarro | 1 | 0.01 |
| Dúplex en Simancas | 2 | 0.03 |
| Dúplex en Trafalgar | 2 | 0.03 |
| Dúplex en Tres Olivos - Valverde | 1 | 0.01 |
| Dúplex en Valdeacederas | 3 | 0.04 |
| Dúplex en Valdebernardo - Valderrivas | 0 | 0.00 |
| Dúplex en Valdemarín | 4 | 0.05 |
| Dúplex en Vallehermoso | 1 | 0.01 |
| Dúplex en Ventilla-Almenara | 0 | 0.00 |
| Dúplex en Villaverde Alto | 1 | 0.01 |
| El Cañaveral | 50 | 0.65 |
| El Pardo | 3 | 0.04 |
| El Plantío | 3 | 0.04 |
| El Viso | 109 | 1.41 |
| Ensanche de Vallecas - La Gavia | 60 | 0.78 |
| Entrevías | 38 | 0.49 |
| Estrella | 27 | 0.35 |
| Fontarrón | 31 | 0.40 |
| Fuente del Berro | 45 | 0.58 |
| Fuentelarreina | 9 | 0.12 |
| Gaztambide | 77 | 1.00 |
| Goya | 257 | 3.32 |
| Guindalera | 97 | 1.25 |
| Hellín | 17 | 0.22 |
| Horcajo | 3 | 0.04 |
| Huertas-Cortes | 48 | 0.62 |
| Ibiza | 56 | 0.72 |
| Imperial | 48 | 0.62 |
| Jerónimos | 62 | 0.80 |
| La Paz | 48 | 0.62 |
| Las Tablas | 39 | 0.50 |
| Lavapiés-Embajadores | 115 | 1.49 |
| Legazpi | 38 | 0.49 |
| Lista | 134 | 1.73 |
| Los Ángeles | 43 | 0.56 |
| Los Berrocales | 3 | 0.04 |
| Los Cármenes | 50 | 0.65 |
| Los Rosales | 34 | 0.44 |
| Lucero | 43 | 0.56 |
| Malasaña-Universidad | 126 | 1.63 |
| Marroquina | 17 | 0.22 |
| Media Legua | 20 | 0.26 |
| Mirasierra | 15 | 0.19 |
| Montecarmelo | 14 | 0.18 |
| Moscardó | 47 | 0.61 |
| Niño Jesús | 24 | 0.31 |
| Nueva España | 76 | 0.98 |
| Nuevos Ministerios-Ríos Rosas | 91 | 1.18 |
| Numancia | 90 | 1.16 |
| Opañel | 54 | 0.70 |
| Orcasitas | 38 | 0.49 |
| Pacífico | 83 | 1.07 |
| Palacio | 88 | 1.14 |
| Palomas | 20 | 0.26 |
| Palomeras Bajas | 42 | 0.54 |
| Palomeras sureste | 37 | 0.48 |
| Palos de Moguer | 73 | 0.94 |
| Pau de Carabanchel | 26 | 0.34 |
| Pavones | 6 | 0.08 |
| Peñagrande | 77 | 1.00 |
| Pilar | 63 | 0.81 |
| Pinar del Rey | 67 | 0.87 |
| Portazgo | 43 | 0.56 |
| Pradolongo | 32 | 0.41 |
| Prosperidad | 86 | 1.11 |
| Pueblo Nuevo | 132 | 1.71 |
| Puerta Bonita | 89 | 1.15 |
| Puerta del Ángel | 137 | 1.77 |
| Quintana | 68 | 0.88 |
| Recoletos | 181 | 2.34 |
| Rejas | 59 | 0.76 |
| Rosas | 30 | 0.39 |
| Salvador | 21 | 0.27 |
| San Cristóbal | 30 | 0.39 |
| San Diego | 120 | 1.55 |
| San Fermín | 33 | 0.43 |
| San Isidro | 91 | 1.18 |
| San Juan Bautista | 39 | 0.50 |
| San Pascual | 48 | 0.62 |
| Sanchinarro | 61 | 0.79 |
| Santa Eugenia | 16 | 0.21 |
| Simancas | 60 | 0.78 |
| Sol | 66 | 0.85 |
| Timón | 14 | 0.18 |
| Trafalgar | 97 | 1.25 |
| Tres Olivos - Valverde | 51 | 0.66 |
| Valdeacederas | 104 | 1.35 |
| Valdebebas - Valdefuentes | 42 | 0.54 |
| Valdebernardo - Valderrivas | 22 | 0.28 |
| Valdemarín | 19 | 0.25 |
| Valdezarza | 42 | 0.54 |
| Vallehermoso | 72 | 0.93 |
| Ventas | 103 | 1.33 |
| Ventilla-Almenara | 46 | 0.59 |
| Villaverde Alto | 87 | 1.13 |
| Vinateros | 20 | 0.26 |
| Virgen del Cortijo - Manoteras | 31 | 0.40 |
| Vista Alegre | 102 | 1.32 |
| Zofío | 25 | 0.32 |
Distrito
En este caso los datos si están correctamente recogidos, repartiendo las casas entre los 21 distritos de Madrid.Casi todos con una representación considerable, en lo que destaca el barrio de Salamanca, con casi el doble de viviendas que el segundo que más tiene, Chamberí.
merge(setNames(as.data.frame(table(mhp_train$district)), c("district", "count")),
setNames(as.data.frame(round(prop.table(table(mhp_train$district))*100, 2)), c("district", "prop (%)"))
) %>%
arrange(desc(count)) %>%
kbl() %>%
kable_material(c("striped", "hover")) %>%
scroll_box(width = "100%", height = "350px")| district | count | prop (%) |
|---|---|---|
| barrio de salamanca | 911 | 11.78 |
| chamberi | 571 | 7.38 |
| ciudad lineal | 540 | 6.98 |
| chamartin | 526 | 6.80 |
| tetuan | 523 | 6.76 |
| centro | 517 | 6.69 |
| carabanchel | 512 | 6.62 |
| puente-de-vallecas | 374 | 4.84 |
| fuencarral | 356 | 4.60 |
| moncloa | 353 | 4.57 |
| hortaleza | 352 | 4.55 |
| latina | 337 | 4.36 |
| san-blas | 319 | 4.13 |
| Arganzuela | 305 | 3.94 |
| retiro | 295 | 3.82 |
| usera | 240 | 3.10 |
| villaverde | 209 | 2.70 |
| villa-de-vallecas | 163 | 2.11 |
| vicalvaro | 141 | 1.82 |
| moratalaz | 99 | 1.28 |
| barajas | 89 | 1.15 |
ggplot(mhp_train, aes(district)) +
geom_bar(fill = "#0BB363") +
labs(x = "Distrito", y = "Número de viviendas", title = "Viviendas por distrito") +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1))
3.2 Análisis de variables cuantitativas
Precio
data.frame(summarise(mhp_train,
min = min(price),
max = max(price),
median = median(price),
mean = mean(price),
sd = sd(price))) %>%
kbl() %>%
kable_material(c("striped", "hover"))| min | max | median | mean | sd |
|---|---|---|---|---|
| 47500 | 8950000 | 369000 | 631920.1 | 753136.6 |
ggplot(mhp_train, aes(y = price)) +
geom_boxplot(fill = "#0BB363") +
labs(y = "Precio", title = "Boxplot de precios") +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1))
ggplot(mhp_train, aes(price)) +
geom_histogram(aes(y = ..count..), bins = 50, position = "dodge", fill = "#0BB363") +
labs(x = "Precio", y = "Número de viviendas", title = "Viviendas por Precio") +
scale_x_continuous(breaks = pretty(mhp_train$price, n = 10), labels = comma_format()) +
scale_y_continuous(labels = function(x) sprintf("%.0f", x)) +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30))
Habitaciones
data.frame(summarise(mhp_train,
min = min(rooms),
max = max(rooms),
median = median(rooms),
mean = mean(rooms),
sd = sd(rooms))) %>%
kbl() %>%
kable_material(c("striped", "hover"))| min | max | median | mean | sd |
|---|---|---|---|---|
| 1 | 17 | 3 | 2.848552 | 1.302995 |
ggplot(mhp_train, aes(y = rooms)) +
geom_boxplot(fill = "#0BB363") +
labs(y = "Habitaciones", title = "Boxplot de habitaciones") +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1))
ggplot(mhp_train, aes(rooms)) +
geom_histogram(aes(), bins = 16, position = "dodge", fill = "#0BB363") +
geom_density(alpha=.2, fill = "red") +
labs(x = "Habitaciones", y = "Número de viviendas", title = "Viviendas por número de habitaciones") +
theme(plot.title = element_text(hjust = 0.5))
m²
data.frame(summarise(mhp_train,
min = min(m2),
max = max(m2),
median = median(m2),
mean = mean(m2),
sd = sd(m2))) %>%
kbl() %>%
kable_material(c("striped", "hover"))| min | max | median | mean | sd |
|---|---|---|---|---|
| 20 | 986 | 95 | 127.1208 | 102.497 |
ggplot(mhp_train, aes(y = m2)) +
geom_boxplot(fill = "#0BB363") +
labs(y = "Metros cuadrados", title = "Boxplot de superficies") +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1))
ggplot(mhp_train, aes(m2)) +
geom_histogram(aes(y=..count..), bins = 50, position = "dodge", fill = "#0BB363") +
geom_density(alpha=.3, fill = "red") +
labs(x = "Metros cuadrados", y = "Número de viviendas", title = "Viviendas por metros cuadrados") +
theme(plot.title = element_text(hjust = 0.5))
4 Análisis multivariante
Una vez analizadas todas las variables de forma individual, podemos buscar algunas relaciones entre ellas.
Como se ve en los gráficos inferiores hay una relción fuerte entre el precio, los metros cuadradados y el número de habitaciones.
plot_grid(
ggcorrplot(cor(mhp_train %>% select_if(is.numeric)),type = "lower", lab=TRUE),
ggplot(mhp_train, aes(x = m2, y = price)) +
geom_point() +
geom_smooth() +
ggtitle('Reación precio y m²') +
theme(plot.title = element_text(hjust = 0.5)),
ggplot(mhp_train, aes(x = rooms, y = price)) +
geom_point() +
geom_smooth() +
ggtitle('Reación precio y habitaciones') +
theme(plot.title = element_text(hjust = 0.5)),
ggplot(mhp_train, aes(x = rooms, y = m2)) +
geom_point() +
geom_smooth() +
ggtitle('Reación m² y habitaciones') +
theme(plot.title = element_text(hjust = 0.5)),
nrow = 2
)
Introduciendo algunas de las variables categóricas se puede observar como varía el precio con la superficie o el número de habitaciones dependiendo de si incluye garaje, ascensor u orientación exterior.
plot_grid(
ggplot(mhp_train, aes(x = m2, y = price, colour = exterior)) +
geom_point() +
geom_smooth(method = "lm") +
scale_colour_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
ggtitle('Rel. precio-superficie, exterior') ,
ggplot(mhp_train, aes(x = rooms, y = price, colour = exterior)) +
geom_point() +
geom_smooth(method = "lm") +
scale_colour_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
ggtitle('Rel. precio-habitaciones, exterior') ,
nrow =1
)
plot_grid(
ggplot(mhp_train, aes(x = m2, y = price, colour = elevator)) +
geom_point() +
geom_smooth(method = "lm") +
scale_colour_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
ggtitle('Rel. precio-superficie, ascensor'),
ggplot(mhp_train, aes(x = rooms, y = price, colour = elevator)) +
geom_point() +
geom_smooth(method = "lm") +
scale_colour_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
ggtitle('Rel. precio-habitación, ascensor'),
nrow =1
)
plot_grid(
ggplot(mhp_train, aes(x = m2, y = price, colour = garage)) +
geom_point() +
geom_smooth(method = "lm") +
scale_colour_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
ggtitle('Rel. precio-superficie, garaje'),
ggplot(mhp_train, aes(x = rooms, y = price, colour = garage)) +
geom_point() +
geom_smooth(method = "lm") +
scale_colour_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
ggtitle('Rel. precio-habitaciones, garaje'),
nrow =1
)
Representando boxplot por distrito se ve claramente que los barrios ricos como Salamanca o Chamberí, a parte de tener precios medios más altos, tienen mayor número de viviendas con precios atípicos (por encima), mientras que en el número de habitaciones no se aprecia una diferencia tan notoria.
ggplot(mhp_train, aes(district, price)) +
geom_boxplot(fill = "#0BB363") +
labs(y = "Precio", title = "Boxplot de precios por distrito") +
scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1))
ggplot(mhp_train, aes(district, rooms)) +
geom_boxplot(fill = "#0BB363") +
labs(y = "Habitaciones", title = "Boxplot de habitaciones por distrito") +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1))
ggplot(mhp_train, aes(district, m2)) +
geom_boxplot(fill = "#0BB363") +
labs(y = "Metros cuadrados", title = "Boxplot de superficie por distrito") +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1))
En las siguientes figuras se representa el precio medio por distrito y como varía con respecto a otras variables de interés.
mhp_train %>%
group_by(district, m2) %>%
summarize(avg_price = mean(price)) %>%
ggplot(aes(x = m2, y = avg_price)) +
geom_point(size = 0.5) +
facet_wrap(~ district) +
scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6))
mhp_train %>%
group_by(district, exterior) %>%
summarise(avg_price = mean(price)) %>%
ggplot(aes(x=district, y=avg_price, fill=exterior)) +
geom_bar(stat = "identity", position = "dodge") +
scale_fill_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
ggtitle("Precio medio por distrito y exterior") +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1)) +
labs(x = "Distrito", y = "Precio medio") +
scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6))
mhp_train %>%
group_by(district, elevator) %>%
summarise(avg_price = mean(price)) %>%
ggplot(aes(x=district, y=avg_price, fill=elevator)) +
geom_bar(stat = "identity", position = "dodge") +
scale_fill_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
ggtitle("Precio medio por distrito y ascensor") +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1)) +
labs(x = "Distrito", y = "Precio medio") +
scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6))
mhp_train %>%
group_by(district, garage) %>%
summarise(avg_price = mean(price)) %>%
ggplot(aes(x=district, y=avg_price, fill=garage)) +
geom_bar(stat = "identity", position = "dodge") +
scale_fill_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
ggtitle("Precio medio por distrito y garaje") +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1)) +
labs(x = "Distrito", y = "Precio medio") +
scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6))
5 Procesado de variables cualitativas
mhp_train %>%
group_by(house_type) %>%
summarise(precio_medio = mean(price)) %>%
arrange(precio_medio)## # A tibble: 26 × 2
## house_type precio_medio
## <fct> <dbl>
## 1 sotano 279031.
## 2 planta -1 298188.
## 3 semi-sotano 302916.
## 4 bajo 316306.
## 5 planta 11 458019.
## 6 planta 1 540375.
## 7 planta 2 595355.
## 8 entreplanta 598834.
## 9 planta 10 605368.
## 10 planta 4 605711.
## # … with 16 more rowsPara el posterior análisis transformamos la columna house_type, agrupando en “piso” (independientemente de la planta en la que esté), “casa” (casa o chalet) y “sotano”.
mhp_train <- mhp_train %>%
mutate(tipo_casa = case_when(
grepl("casa|chalet", house_type, ignore.case = TRUE) ~ "casa",
grepl("\\bplanta\\s*-?1\\b|semi-?sotano|bajo|sotano", house_type, ignore.case = TRUE) ~ "sotano",
TRUE ~ "piso"
))
mhp_train$tipo_casa <- as.factor(mhp_train$tipo_casa)También se agruparán los distritos en función del precio medio del mismo: caro, medio o barato.
mhp_train <- mhp_train %>%
group_by(district) %>%
mutate(precio_distrito = mean(price)) %>%
ungroup()
# Calcular los percentiles
percentiles <- quantile(mhp_train$precio_distrito, probs = c(0, 1/3, 2/3, 1))
# Asignar etiquetas
mhp_train$precio_distrito <- cut(mhp_train$precio_distrito, breaks = percentiles, labels = c("barato", "medio", "caro"), include.lowest = TRUE)Además, de cara a desarrollar algunos modelos, se creará una variable target binaria que dividirá las casas entre “cara” y “barata” (1/0).
mhp_train$precio_bin <- ifelse(mhp_train$price > median(mhp_train$price), "cara", "barata")
mhp_train$precio_bin = ifelse(mhp_train$precio_bin == "cara", 1, 0)
mhp_train$precio_bin = factor(mhp_train$precio_bin, levels = c(1,0))# Aplicamos todos los cambios también a test.
# Categorizamos house_type
mhp_test <- mhp_test %>%
mutate(tipo_casa = case_when(
grepl("casa|chalet", house_type, ignore.case = TRUE) ~ "casa",
grepl("\\bplanta\\s*-?1\\b|semi-?sotano|bajo|sotano", house_type, ignore.case = TRUE) ~ "sotano",
TRUE ~ "piso"
))
mhp_test$tipo_casa <- as.factor(mhp_test$tipo_casa)
# Categorizamos los distritos
mhp_test <- mhp_test %>%
group_by(district) %>%
mutate(precio_distrito = mean(price)) %>%
ungroup()
percentiles <- quantile(mhp_test$precio_distrito, probs = c(0, 1/3, 2/3, 1))
mhp_test$precio_distrito <- cut(mhp_test$precio_distrito, breaks = percentiles, labels = c("barato", "medio", "caro"), include.lowest = TRUE)
# Creamos variable target
mhp_test$precio_bin <- ifelse(mhp_test$price > median(mhp_test$price), "cara", "barata")
mhp_test$precio_bin = ifelse(mhp_test$precio_bin == "cara", 1, 0)
mhp_test$precio_bin = factor(mhp_test$precio_bin, levels = c(1,0))6 Modelo de regresión lineal
Con el dataset ya preparado, se crea un modelo de regresión lineal.
En este caso, tras hacer diversas pruebas se ha optado por introducir sólo 3 variables (tipo_casa, m2 y precio_distrito), alcanzando un R2 de 0.7553.
Esta cifra podía incrementarse algo, sin embargo es a costa de introducir más variables y por tanto más complejidad al modelo, haciéndolo menos explicable.
lm_fit <- lm(price ~ m2+tipo_casa+precio_distrito, data=mhp_train)
summary(lm_fit)##
## Call:
## lm(formula = price ~ m2 + tipo_casa + precio_distrito, data = mhp_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3277290 -137875 -17495 95705 5483233
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -915949.86 30574.68 -29.958 <2e-16 ***
## m2 6052.79 51.28 118.027 <2e-16 ***
## tipo_casapiso 673879.37 27865.36 24.183 <2e-16 ***
## tipo_casasotano 606283.79 28773.95 21.071 <2e-16 ***
## precio_distritomedio 101183.27 10337.15 9.788 <2e-16 ***
## precio_distritocaro 373747.30 11419.31 32.729 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 372600 on 7726 degrees of freedom
## Multiple R-squared: 0.7553, Adjusted R-squared: 0.7552
## F-statistic: 4771 on 5 and 7726 DF, p-value: < 2.2e-16coef(lm_fit)## (Intercept) m2 tipo_casapiso
## -915949.864 6052.795 673879.371
## tipo_casasotano precio_distritomedio precio_distritocaro
## 606283.787 101183.273 373747.298residuals=lm_fit$residuals
autoplot(lm_fit)
Observando los residuos se ve como hay cierta heterocedasticidad y en la gráfica QQ no se ajustan a la normal, por lo que el módelo no es bueno.
# Realizar predicciones en la partición test
predictions <- predict(lm_fit, newdata = mhp_test)
# Comparar las predicciones con los valores reales
rmse <- rmse(predictions, mhp_test$price)
mae <- mae(predictions, mhp_test$price)
r2 <- R2(predictions, mhp_test$price)
# Imprimir las métricas de evaluación
cat(paste0("RMSE: ", round(rmse, 2), "\n"))## RMSE: 398263.79cat(paste0("MAE: ", round(mae, 2), "\n"))## MAE: 199633.64cat(paste0("R-squared: ", round(r2, 2)))## R-squared: 0.73Un R2 de 0.73 (bastante similar al calculado en train) indica que se explica con esa probabilidad la variabilidad de la variable dependiente, por lo tanto podríamos considerarlo razonablemente útil, sin embargo sería preferible mejorar los valores del error de predicción promedio o el cuadrático.
7 Conclusiones preliminares
El análisis realizado hasta el momento ha servido para familiarizarse con el mercado inmobiliario de Madrid (o al menos parte de él) y conocer como se comportan algunas de las variables más relevantes en el valor de las viviendas.
Con estos conocimientos se abordarán nuevos modelos para intentar conseguir mejores predicciones y simultaneamente comprender y desarrollar estos mismos modelos.
str(mhp_train)## tibble [7,732 × 12] (S3: tbl_df/tbl/data.frame)
## $ price : num [1:7732] 499000 385000 950000 460000 335000 580000 940000 85300 399000 125000 ...
## $ house_type : Factor w/ 26 levels "bajo","casa",..: 21 16 20 10 18 19 16 16 6 1 ...
## $ exterior : Factor w/ 2 levels "1","0": 1 1 1 1 1 1 1 1 1 1 ...
## $ rooms : num [1:7732] 4 1 3 3 4 3 2 3 3 2 ...
## $ m2 : num [1:7732] 131 60 160 114 122 120 124 75 86 62 ...
## $ elevator : Factor w/ 2 levels "1","0": 1 2 1 1 1 1 1 1 2 2 ...
## $ garage : Factor w/ 2 levels "1","0": 2 2 2 1 2 2 2 2 2 2 ...
## $ neighborhood : Factor w/ 341 levels "12 de Octubre-Orcasur",..: 76 17 8 284 318 265 264 316 192 11 ...
## $ district : Factor w/ 21 levels "Arganzuela","barajas",..: 17 12 7 13 18 3 3 21 1 19 ...
## $ tipo_casa : Factor w/ 3 levels "casa","piso",..: 2 2 2 2 2 2 2 2 3 3 ...
## $ precio_distrito: Factor w/ 3 levels "barato","medio",..: 2 3 3 1 1 3 3 1 1 1 ...
## $ precio_bin : Factor w/ 2 levels "1","0": 1 1 1 1 2 1 1 2 1 2 ...8 Aprendizaje no supervisado
8.1 K-Means
El algoritmo K-Means es un método de aprendizaje automático no supervisado utilizado para el análisis de agrupamiento o “clustering”. El objetivo principal de este algoritmo es dividir los datos en K grupos (clusters), de modo que los puntos dentro de un mismo grupo sean lo más similares posible.
Por como funciona este método, es habitual utilizarlo en la etapa de EDA para agrupar y etiquetar observaciones con un criterior más científico que el basado en el conocimiento personal del dominio, más que para hacer predicciones (aunque también puede hacerse).
En nuestro lo utilizaremos para comprobar si hay alguna forma mejor de agrupar las casas de como lo hemos hecho anteriormente.
mhp_train_kmeans <- subset(mhp_train, select = c("m2", "rooms", "elevator", "garage", "exterior"))
mhp_train_kmeans$exterior <- as.numeric(mhp_train_kmeans$exterior)
mhp_train_kmeans$elevator <- as.numeric(mhp_train_kmeans$elevator)
mhp_train_kmeans$garage <- as.numeric(mhp_train_kmeans$garage)
# Normalizamos los datos
mhp_train_kmeans <- scale(mhp_train_kmeans)Tras normalizar los datos, aplicamos k-means sólo para las variables numéricas y mediante los siguientes test decidiremos cual es el número óptimo de clusters.
fviz_nbclust(mhp_train_kmeans, kmeans, method = "wss")
fviz_nbclust(mhp_train_kmeans, kmeans, method = "silhouette")
Utilizando la regla del codo en el gráfico wss podríamos decir que está en el 6, algo que confirma el segundo gráfico. Por lo tanto, aplicamos el modelo con K = 6
# Definir el número de grupos
k <- 6
# Ejecutar k-means
grupos <- kmeans(mhp_train_kmeans, k)
mhp_train$cluster <- as.factor(grupos$cluster)
# Graficar la relación entre las características y los grupos
ggplot(mhp_train, aes(x=m2, y=price, color=cluster)) +
geom_point() +
ggtitle("Grupos de casas en venta") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Tamaño en m2") +
ylab("Precio")
fviz_cluster(grupos, data = mhp_train_kmeans, ellipse.type = "euclid",repel = TRUE,star.plot = TRUE)
fviz_cluster(grupos, data = mhp_train_kmeans, ellipse.type = "norm")
Graficamente es difícil diferenciar las agrupaciones, por lo que en la tabla inferior obtenemos los precios medios de cada cluster, donde se puede apreciar que la gran mayoría de agrupaciones son coherentes (salvo tal vez el 2 y el 3 que son muy similares).
También se muestra el tipo de casa y el distrito predominante en cada cluster.
# Precio medio por cluster
mhp_train %>%
group_by(cluster) %>%
summarise(precio_medio = mean(price))## # A tibble: 6 × 2
## cluster precio_medio
## <fct> <dbl>
## 1 1 184203.
## 2 2 640443.
## 3 3 617779.
## 4 4 223476.
## 5 5 2228444.
## 6 6 355840.# Tipo de casa predominante por cluster
mhp_train %>%
group_by(cluster, tipo_casa) %>%
summarise(frecuencia = n()) %>%
arrange(cluster, desc(frecuencia)) %>%
slice(1) %>%
select(cluster, tipo_casa_predominante = tipo_casa)## # A tibble: 6 × 2
## # Groups: cluster [6]
## cluster tipo_casa_predominante
## <fct> <fct>
## 1 1 sotano
## 2 2 piso
## 3 3 piso
## 4 4 piso
## 5 5 piso
## 6 6 piso# Distrito predominante por cluster
mhp_train %>%
group_by(cluster, district) %>%
summarise(frecuencia = n()) %>%
arrange(cluster, desc(frecuencia)) %>%
slice(1) %>%
select(cluster, distrito_predominante = district)## # A tibble: 6 × 2
## # Groups: cluster [6]
## cluster distrito_predominante
## <fct> <fct>
## 1 1 tetuan
## 2 2 hortaleza
## 3 3 barrio de salamanca
## 4 4 carabanchel
## 5 5 barrio de salamanca
## 6 6 barrio de salamancaSabemos que tanto la regla del codo como cualquier otro método para definir el número de clusters no siempre es adecuada, por lo que a continuación se prueba con un bucle que prueba con número de clusters menores que 6 y muestra las agrupaciones por precio y tipo de piso anteriores.
# Ejecuta k-means para diferentes valores de K
for (k in 2:6) {
# Ejecutar k-means
grupos <- kmeans(mhp_train_kmeans, centers = k)
mhp_train$cluster <- as.factor(grupos$cluster)
print(k)
print(mhp_train %>%
group_by(cluster) %>%
summarise(precio_medio = mean(price)))
print(mhp_train %>%
group_by(cluster, tipo_casa) %>%
summarise(frecuencia = n()) %>%
arrange(cluster, desc(frecuencia)) %>%
slice(1) %>%
select(cluster, tipo_casa_predominante = tipo_casa))
}## [1] 2
## # A tibble: 2 × 2
## cluster precio_medio
## <fct> <dbl>
## 1 1 394442.
## 2 2 1094473.
## # A tibble: 2 × 2
## # Groups: cluster [2]
## cluster tipo_casa_predominante
## <fct> <fct>
## 1 1 piso
## 2 2 piso
## [1] 3
## # A tibble: 3 × 2
## cluster precio_medio
## <fct> <dbl>
## 1 1 402986.
## 2 2 605000.
## 3 3 1941976.
## # A tibble: 3 × 2
## # Groups: cluster [3]
## cluster tipo_casa_predominante
## <fct> <fct>
## 1 1 piso
## 2 2 piso
## 3 3 piso
## [1] 4
## # A tibble: 4 × 2
## cluster precio_medio
## <fct> <dbl>
## 1 1 579636.
## 2 2 352563.
## 3 3 665255.
## 4 4 189057.
## # A tibble: 4 × 2
## # Groups: cluster [4]
## cluster tipo_casa_predominante
## <fct> <fct>
## 1 1 piso
## 2 2 piso
## 3 3 piso
## 4 4 sotano
## [1] 5
## # A tibble: 5 × 2
## cluster precio_medio
## <fct> <dbl>
## 1 1 435798.
## 2 2 241929.
## 3 3 189057.
## 4 4 665853.
## 5 5 1005091.
## # A tibble: 5 × 2
## # Groups: cluster [5]
## cluster tipo_casa_predominante
## <fct> <fct>
## 1 1 piso
## 2 2 piso
## 3 3 sotano
## 4 4 piso
## 5 5 piso
## [1] 6
## # A tibble: 6 × 2
## cluster precio_medio
## <fct> <dbl>
## 1 1 1179261.
## 2 2 652498.
## 3 3 355840.
## 4 4 217255.
## 5 5 2356279.
## 6 6 451222.
## # A tibble: 6 × 2
## # Groups: cluster [6]
## cluster tipo_casa_predominante
## <fct> <fct>
## 1 1 piso
## 2 2 piso
## 3 3 piso
## 4 4 piso
## 5 5 piso
## 6 6 pisoObservando las diferencias entre los diferentes modelos, el número de clusters que mejor diferencia las casas por su precio es 3, algo que encaja con la presunción que inicial que hicimos para dividir los tipos de casa en 3 tipos.
Si observamos además el tipo de casas más común por cluster, casi siempre es piso (salvo en los casos de clusters con muy poco precio medio, que es sótano), lo cual aunque tiene sentido porque este tipo es claramente mayoritario en todo el dataset, no nos permite apreciar a simple vista donde se acumulan el resto de viviendas.
Por lo tanto, entre los datos extraídos del análisis de clusters y el conocimiento previo que teníamos del dominio, tiene sentido mantener los 3 tipos de viviendas principales como variable para los siguientes modelos (sótano, piso y casa).
9 Técnicas de reducción de la dimensionalidad
9.1 PCA
El Análisis de Componentes Principales es una técnica de reducción de la dimensionalidad que se utiliza para transformar un conjunto de datos de múltiples variables en un nuevo conjunto de datos con un número menor de variables cuyo principal objetivo es simplificar la representación de datos mientras se conserva la mayor cantidad de información y variabilidad posible identificando patrones y estructuras subyacentes, reduciendo la multicolinealidad y el sobreajuste.
PAra nuestro dataset, seleccionaremos las columnas numéricas, descartando la columna precio puesto que será de alguna manera nuestra variable dependiente. Además creamos 3 variables dummy usando variables categóricas de tipo 0|1.
train_pca <- mhp_train[,c(3:7,10,11,12)]dummy_garaje <- model.matrix(~0 + garage,data = train_pca)
dummy_exterior <- model.matrix(~0 + exterior, data = train_pca)
dummy_ascensor <- model.matrix(~0 + elevator, data = train_pca)
dummy_tipo_casa <- model.matrix(~0 + tipo_casa, data = train_pca)
dummy_distrito <- model.matrix(~0 + precio_distrito, data = train_pca)dummy_data <- cbind(train_pca[, c("rooms", "m2")],dummy_exterior,dummy_garaje,dummy_ascensor,dummy_distrito,dummy_tipo_casa)prcomp(dummy_data)## Standard deviations (1, .., p=14):
## [1] 1.025015e+02 9.457105e-01 7.267966e-01 6.257613e-01 6.127629e-01
## [6] 5.782676e-01 4.571490e-01 4.037562e-01 1.716715e-01 2.268386e-15
## [11] 1.853843e-15 1.003401e-15 8.037840e-16 4.438592e-16
##
## Rotation (n x k) = (14 x 14):
## PC1 PC2 PC3 PC4
## rooms -8.797027e-03 -0.9916184565 -0.014616928 0.094846395
## m2 -9.999557e-01 0.0085416356 0.001971791 -0.001450128
## exterior1 -5.771447e-04 -0.0453751425 -0.009449951 -0.063478191
## exterior0 5.771447e-04 0.0453751425 0.009449951 0.063478191
## garage1 -1.538147e-03 0.0421269404 -0.223561300 0.037355645
## garage0 1.538147e-03 -0.0421269404 0.223561300 -0.037355645
## elevator1 -4.274374e-04 0.0241411854 -0.476833313 0.342257253
## elevator0 4.274374e-04 -0.0241411854 0.476833313 -0.342257253
## precio_distritobarato 1.400700e-03 -0.0385563529 0.200000119 -0.359200946
## precio_distritomedio 6.834450e-05 -0.0051663464 -0.051807433 0.201186914
## precio_distritocaro -1.469044e-03 0.0437226993 -0.148192687 0.158014031
## tipo_casacasa -8.085078e-04 0.0005892036 0.067961955 -0.037840260
## tipo_casapiso 2.818685e-05 -0.0469938502 -0.466029098 -0.505128722
## tipo_casasotano 7.803209e-04 0.0464046466 0.398067143 0.542968982
## PC5 PC6 PC7 PC8
## rooms -0.0132115607 -0.0291253564 0.074126853 -0.029792435
## m2 0.0009596296 0.0005665667 -0.002313108 -0.001386419
## exterior1 -0.1552630611 -0.0333936354 -0.327138796 0.601052285
## exterior0 0.1552630611 0.0333936354 0.327138796 -0.601052285
## garage1 -0.5627134868 -0.1309300044 0.333444368 0.033280690
## garage0 0.5627134868 0.1309300044 -0.333444368 -0.033280690
## elevator1 0.0482119831 -0.1533456743 -0.318624110 -0.141577832
## elevator0 -0.0482119831 0.1533456743 0.318624110 0.141577832
## precio_distritobarato -0.2366848900 -0.4502445095 -0.365221920 -0.322493702
## precio_distritomedio -0.1969164022 0.7585991702 -0.088292084 -0.036423761
## precio_distritocaro 0.4336012922 -0.3083546607 0.453514004 0.358917463
## tipo_casacasa -0.0381883034 0.0229829599 0.061657352 -0.001868362
## tipo_casapiso 0.1388358808 0.1326235059 0.004437444 -0.006913037
## tipo_casasotano -0.1006475774 -0.1556064657 -0.066094796 0.008781399
## PC9 PC10 PC11 PC12
## rooms 0.000862928 5.877890e-16 1.174221e-15 7.345674e-17
## m2 -0.001036478 -4.349751e-20 3.574359e-18 -4.762673e-18
## exterior1 0.017792902 8.077407e-03 2.416838e-02 -9.348594e-03
## exterior0 -0.017792902 8.077407e-03 2.416838e-02 -9.348594e-03
## garage1 -0.027676182 6.464949e-03 5.102210e-03 -8.893342e-02
## garage0 0.027676182 6.464949e-03 5.102210e-03 -8.893342e-02
## elevator1 0.086614598 6.353245e-01 -3.028358e-01 -6.708710e-02
## elevator0 -0.086614598 6.353245e-01 -3.028358e-01 -6.708710e-02
## precio_distritobarato -0.004863213 -9.943406e-02 -8.213108e-02 -5.584679e-01
## precio_distritomedio -0.010432322 -9.943406e-02 -8.213108e-02 -5.584679e-01
## precio_distritocaro 0.015295535 -9.943406e-02 -8.213108e-02 -5.584679e-01
## tipo_casacasa 0.809211511 -2.329879e-01 -5.148217e-01 1.145226e-01
## tipo_casapiso -0.393939361 -2.329879e-01 -5.148217e-01 1.145226e-01
## tipo_casasotano -0.415272150 -2.329879e-01 -5.148217e-01 1.145226e-01
## PC13 PC14
## rooms -9.651041e-17 -5.224615e-18
## m2 4.139246e-18 -4.009145e-19
## exterior1 5.278108e-02 -7.046115e-01
## exterior0 5.278108e-02 -7.046115e-01
## garage1 -6.995888e-01 -5.097578e-02
## garage0 -6.995888e-01 -5.097578e-02
## elevator1 1.228500e-02 -1.293890e-03
## elevator0 1.228500e-02 -1.293890e-03
## precio_distritobarato 6.884856e-02 8.609917e-03
## precio_distritomedio 6.884856e-02 8.609917e-03
## precio_distritocaro 6.884856e-02 8.609917e-03
## tipo_casacasa -1.877163e-02 -2.325502e-02
## tipo_casapiso -1.877163e-02 -2.325502e-02
## tipo_casasotano -1.877163e-02 -2.325502e-02Podemos observar que en la PC1 lleva un peso casi de 100 la variable m2, mientras que en PC2 pasa lo mismo con la variable rooms. En el resto de PCs el análisis de componentes principales está algo más repartido.
summary(prcomp(dummy_data))## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 102.5015 0.94571 0.72680 0.62576 0.61276 0.57827 0.45715
## Proportion of Variance 0.9997 0.00009 0.00005 0.00004 0.00004 0.00003 0.00002
## Cumulative Proportion 0.9997 0.99981 0.99986 0.99989 0.99993 0.99996 0.99998
## PC8 PC9 PC10 PC11 PC12 PC13
## Standard deviation 0.40376 0.1717 2.268e-15 1.854e-15 1.003e-15 8.038e-16
## Proportion of Variance 0.00002 0.0000 0.000e+00 0.000e+00 0.000e+00 0.000e+00
## Cumulative Proportion 1.00000 1.0000 1.000e+00 1.000e+00 1.000e+00 1.000e+00
## PC14
## Standard deviation 4.439e-16
## Proportion of Variance 0.000e+00
## Cumulative Proportion 1.000e+00Al hacer un summary de las variables vemos que la proporción de varianza explicada por la primera componente es de un 99%, es decir, que es la única relevante. Si nos fijamos vemos que la PC1 está explicada casi al 100% por los m2. En nuestra tabla vemos que los m2 tienen valores mucho más altos que el resto de variables. Por lo tanto, lo que hacemos a continuación es un análisis pero escalando las variables.
prcomp(dummy_data, scale = T)## Standard deviations (1, .., p=14):
## [1] 1.788828e+00 1.554767e+00 1.391248e+00 1.320487e+00 1.255039e+00
## [6] 1.218340e+00 9.006738e-01 7.679348e-01 4.930742e-01 2.297198e-14
## [11] 4.169026e-15 2.788355e-15 1.325148e-15 8.186691e-16
##
## Rotation (n x k) = (14 x 14):
## PC1 PC2 PC3 PC4
## rooms -0.32118124 0.185174155 0.19374853 -0.18231760
## m2 -0.39188093 0.178866531 0.28992752 -0.15926158
## exterior1 -0.22525556 0.323113180 -0.41045579 0.16000993
## exterior0 0.22525556 -0.323113180 0.41045579 -0.16000993
## garage1 -0.39123327 0.104563202 0.02079773 0.20853492
## garage0 0.39123327 -0.104563202 -0.02079773 -0.20853492
## elevator1 -0.28418360 -0.450986811 -0.07378472 0.27110372
## elevator0 0.28418360 0.450986811 0.07378472 -0.27110372
## precio_distritobarato 0.20341246 0.208545636 -0.33198104 0.01261772
## precio_distritomedio -0.02774634 -0.002015776 0.03875669 0.15703618
## precio_distritocaro -0.18188148 -0.213762808 0.30358073 -0.17520956
## tipo_casacasa -0.14205765 0.358486679 0.29646454 -0.20249803
## tipo_casapiso -0.16858177 -0.251362769 -0.39199358 -0.48519863
## tipo_casasotano 0.22431380 0.124024528 0.29030539 0.56945631
## PC5 PC6 PC7 PC8
## rooms 0.17534372 -0.091653082 0.472849806 0.496683083
## m2 0.12920863 -0.038990789 0.209138336 0.106709918
## exterior1 0.31433452 -0.162497819 -0.145528928 -0.064746601
## exterior0 -0.31433452 0.162497819 0.145528928 0.064746601
## garage1 -0.40131726 0.276379019 -0.203710793 0.095975695
## garage0 0.40131726 -0.276379019 0.203710793 -0.095975695
## elevator1 0.11627307 0.012675714 0.278402861 -0.212854201
## elevator0 -0.11627307 -0.012675714 -0.278402861 0.212854201
## precio_distritobarato -0.04064647 0.512859981 0.434958282 -0.014241486
## precio_distritomedio -0.38476634 -0.695032433 0.035559829 0.020880654
## precio_distritocaro 0.43936573 0.186858173 -0.486898651 -0.006820754
## tipo_casacasa -0.08089606 -0.001818982 0.145439408 -0.763400011
## tipo_casapiso -0.15088701 -0.039033227 -0.059791130 0.108495318
## tipo_casasotano 0.18370301 0.040471673 0.007300261 0.171061449
## PC9 PC10 PC11 PC12
## rooms -0.53136448 1.871326e-14 6.573067e-16 -8.713797e-16
## m2 0.79477851 5.562477e-15 5.135097e-16 5.962879e-16
## exterior1 -0.01407624 -1.579692e-03 -7.045930e-01 -5.852591e-02
## exterior0 0.01407624 -1.579692e-03 -7.045930e-01 -5.852591e-02
## garage1 -0.06271669 -7.002484e-04 -1.864879e-03 2.225429e-02
## garage0 0.06271669 -7.002484e-04 -1.864879e-03 2.225429e-02
## elevator1 -0.02034701 -1.325824e-03 -1.363473e-02 3.145845e-02
## elevator0 0.02034701 -1.325824e-03 -1.363473e-02 3.145845e-02
## precio_distritobarato 0.09621994 -7.055211e-02 4.768696e-02 -5.771886e-01
## precio_distritomedio 0.02039646 -7.038638e-02 4.757494e-02 -5.758327e-01
## precio_distritocaro -0.12064810 -6.816645e-02 4.607446e-02 -5.576714e-01
## tipo_casacasa -0.22215925 -2.483203e-01 -1.912383e-03 3.007261e-02
## tipo_casapiso 0.04299775 -6.861954e-01 -5.284580e-03 8.310109e-02
## tipo_casasotano 0.03813232 -6.729672e-01 -5.182706e-03 8.149911e-02
## PC13 PC14
## rooms -4.757936e-17 3.988302e-17
## m2 -1.440589e-18 -3.087024e-16
## exterior1 1.099278e-02 3.536496e-04
## exterior0 1.099278e-02 3.536496e-04
## garage1 2.157713e-02 -7.064242e-01
## garage0 2.157713e-02 -7.064242e-01
## elevator1 -7.059752e-01 -2.053508e-02
## elevator0 -7.059752e-01 -2.053508e-02
## precio_distritobarato -2.595459e-02 -1.903173e-02
## precio_distritomedio -2.589362e-02 -1.898702e-02
## precio_distritocaro -2.507696e-02 -1.838819e-02
## tipo_casacasa 1.806856e-03 1.253756e-03
## tipo_casapiso 4.992972e-03 3.464565e-03
## tipo_casasotano 4.896719e-03 3.397776e-03Tras escalar las variables podemos observar que están mucho más repartidos los pesos en las componentes principales debido a que estás asegurando que todas las variables tengan la misma importancia en el análisis, independientemente de sus escalas originales y unidades de medida. Con esto evitas el dominio de variables con mayor magnitud y mejoras la interpretación de resultados.
pca_result <- prcomp(dummy_data, center = TRUE, scale = TRUE)
summary(pca_result)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.7888 1.5548 1.3912 1.3205 1.2550 1.2183 0.90067
## Proportion of Variance 0.2286 0.1727 0.1383 0.1245 0.1125 0.1060 0.05794
## Cumulative Proportion 0.2286 0.4012 0.5395 0.6640 0.7765 0.8826 0.94051
## PC8 PC9 PC10 PC11 PC12 PC13
## Standard deviation 0.76793 0.49307 2.297e-14 4.169e-15 2.788e-15 1.325e-15
## Proportion of Variance 0.04212 0.01737 0.000e+00 0.000e+00 0.000e+00 0.000e+00
## Cumulative Proportion 0.98263 1.00000 1.000e+00 1.000e+00 1.000e+00 1.000e+00
## PC14
## Standard deviation 8.187e-16
## Proportion of Variance 0.000e+00
## Cumulative Proportion 1.000e+00Observamos que con PC6 ya llegamos al 88% de la proporción de varianza explicada pero vemos algunos saltos interesantes. Para ver cuantas componentes explican el modelo hacemos la regla de codo, método heurístico para determinar el número óptimo de clusters en un conjunto de datos.
Al observar el gráfico vemos que la forma del codo sería en la 6|7 aunque el gráfico no es muy claro.
summary_pca <- summary(pca_result)
var_exp <- summary_pca$importance[2,]
barplot(var_exp, col = "#0BB363",
main = "Gráfico PCA - Variance explained",
xlab = "Principal Components",
ylab = "Proportion of Variance Explained")
lines(var_exp, col = "#DC143C", type = "l", lwd = 2)
points(var_exp, col = "#DC143C", pch = 20, cex = 1.2)
A modo de ejercicio para comprender el funcionamiento de PCA puede resultar útil hacer estos cálculos, sin embargo para el dataset en cuestión ya se pueden formar modelos robustos con muy pocas variables (los m2 de superficie, por ejemplo, ya son suficientes para predecir un gran porcentaje de precios).
10 Aprendizaje supervisado
10.1 GLM
En aprendizaje automático, GLM (Generalized Linear Models) se refiere a una clase de modelos lineales que se utilizan para predecir una variable respuesta a partir de una o más variables explicativas. Tiene tres componentes: Variable respuesta, función de enlace y función lineal sistemática, y son muy útiles para utilizarse en clasificación, regresión y problemas de predicción.
Aplicado a nuestro dataset, lo utilizaremos para predecir el precio binario de las casas (cara/barata).
train_glm <- mhp_train[,c(3:7,10,11,12)]
test_glm <- mhp_test[,c(3:7,10,11,12)]train_glm[,2:3]<- scale(train_glm[,2:3])
test_glm[,2:3]<- scale(test_glm[,2:3])Tras seleccionar las variables y normalizarlas para que todas tengan el mismo peso, desarrollamos el modelo con la variable de respuesta “precio_bin” y las variables independientes “exterior”, “garage”, “elevator”, “rooms”, “precio_distrito” y “tipo_casa”, usando una distribución binomial.
glm_mhp = glm(formula = precio_bin ~ exterior + garage + elevator + rooms + m2 + precio_distrito + tipo_casa,
data = train_glm,
family = binomial)
summary(glm_mhp)##
## Call:
## glm(formula = precio_bin ~ exterior + garage + elevator + rooms +
## m2 + precio_distrito + tipo_casa, family = binomial, data = train_glm)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.5396 -0.1705 0.0219 0.2836 5.0687
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.99261 0.62333 -1.592 0.11129
## exterior0 -0.35269 0.13261 -2.660 0.00782 **
## garage0 0.79113 0.10737 7.368 1.73e-13 ***
## elevator0 1.67266 0.12813 13.054 < 2e-16 ***
## rooms -0.01548 0.08151 -0.190 0.84940
## m2 -6.20689 0.23649 -26.246 < 2e-16 ***
## precio_distritomedio -2.52394 0.11266 -22.403 < 2e-16 ***
## precio_distritocaro -4.33905 0.14476 -29.974 < 2e-16 ***
## tipo_casapiso 0.34751 0.62328 0.558 0.57716
## tipo_casasotano 1.01944 0.62491 1.631 0.10282
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 10718.8 on 7731 degrees of freedom
## Residual deviance: 3535.9 on 7722 degrees of freedom
## AIC: 3555.9
##
## Number of Fisher Scoring iterations: 8En este caso, todos los coeficientes tienen valores “p” muy pequeños, lo que indica que todas las variables independientes tienen una relación significativa con la variable dependiente “precio_bin” en el modelo de regresión logística ajustado.
Siendo el “Intercept” negativo, se puede apreciar como el precio depende de forma proporcional con el número de habitaciones, el tipo de piso y los metros cuadrados (sobre todo). Todo esperable, así como la dependencia inversa con aquellas viviendas que no tienen ascensor o graje.
El modelo explica una cantidad significativa de la varianza en la variable respuesta, con un AIC de 4818.
head(predict(glm_mhp))## 1 2 3 4 5 6
## -2.6264925 1.5662145 -6.1858748 0.1476553 0.4424557 -3.7636044Esto sirve para predecir la probabilidad logit ajustado previamente y para mostrar las primeras seis predicciones. Cada número representa la probabilidad logit de que la variable dependiente “precio_bin” sea igual a 1 para cada observación en el conjunto de datos.
ajustados <- fitted(glm_mhp)
head(fitted(glm_mhp))## 1 2 3 4 5 6
## 0.067452746 0.827243283 0.002054072 0.536846908 0.608844012 0.022673932Aquí obtenemos las probabilidades ajustadas que representan las probabilidades estimadas de que la variable dependiente “precio_bin” sea igual a 1 para cada observación en el conjunto de datos.
Gráficos de residuales: Estos gráficos permiten evaluar la calidad del ajuste del modelo. Puedes trazar residuales de Pearson, deviance, o residuales estandarizados frente a las predicciones ajustadas o valores ajustados para examinar si hay patrones en los residuales. Un buen ajuste del modelo no mostrará patrones claros en los residuales.
residuales <- residuals(glm_mhp, type = "pearson")
head(residuales)## 1 2 3 4 5 6
## -0.26894557 -2.18826117 -0.04536849 -1.07662124 0.80153404 -0.15231535data_residuales <- data.frame(Residuales = residuales, Ajustados = ajustados)
ggplot(data_residuales, aes(x = Ajustados, y = Residuales)) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(title = "Gráfico de residuales", x = "Valores ajustados", y = "Residuales") + ylim(-10,10)
predicciones_glm <- ifelse(test = glm_mhp$fitted.values > 0.5, no = 0, yes = 1)
matriz_confusion_glm <- table(glm_mhp$model$precio_bin, predicciones_glm,
dnn = c("observaciones", "predicciones"))
matriz_confusion_glm## predicciones
## observaciones 0 1
## 1 3449 414
## 0 329 3540Finalmente generamos predicciones binarias basadas en las probabilidades ajustadas del modelo lineal generalizado y creamos una matriz de confusión para evaluar el rendimiento del modelo donde se muestra la distribución de las predicciones correctas e incorrectas en función de las observaciones reales y las predichas.
La matriz se lee de la siguiente manera:
Verdaderos negativos (VN): La esquina superior izquierda (3265) representa el número de observaciones clasificadas correctamente como “0” (clase negativa). Esto significa que el modelo predijo correctamente 3265 casos como “0”.
Falsos positivos (FP): La esquina superior derecha (598) representa el número de observaciones que en realidad eran “0” pero fueron clasificadas incorrectamente como “1” (clase positiva) por el modelo. El modelo cometió un error en 598 casos al clasificarlos como “1”.
Falsos negativos (FN): La esquina inferior izquierda (501) representa el número de observaciones que en realidad eran “1” pero fueron clasificadas incorrectamente como “0” por el modelo. El modelo cometió un error en 501 casos al clasificarlos como “0”.
Verdaderos positivos (VP): La esquina inferior derecha (3368) representa el número de observaciones clasificadas correctamente como “1”. Esto significa que el modelo predijo correctamente 3368 casos como “1”.
10.2 KNN
El modelo KNN (K-Nearest Neighbors) es un algoritmo de aprendizaje supervisado que se utiliza tanto para clasificación como para regresión. Es un método basado en instancias que no hace supuestos sobre la forma funcional de la relación entre las variables independientes y la variable dependiente. KNN utiliza la información de los vecinos más cercanos para realizar predicciones. Para ello calcula distancias como la Euclidiana, Manhattan, Minkowski, etc…
La idea detrás de KNN es simple, una observación se clasifica o se le asigna un valor basado en las características de sus K vecinos más cercanos.
Algunas ventajas del algoritmo KNN son su simplicidad o su capacidad para manejar relaciones no lineales entre variables. Sin embargo también tiene algunas desventajas, como su sensibilidad a la maldición de la dimensionalidad, el efecto de las variables irrelevantes y el costo computacional asociado con el cálculo de distancias en conjuntos de datos grandes.
Procedemos al cálculo de este modelo e intentamos ver si hay un número correcto de k-vecinos que podemos usar.
long = 15
accuracy = rep(0,long)
f1score = rep(0,long)
recall = rep(0,long)
precision = rep(0,long)
for (i in 1:long)
{
prediccion_knn_cv =knn.cv(mhp_train[,c("exterior","rooms","m2","elevator", "garage")],
k=i, cl=mhp_train$precio_bin)
accuracy[i] = sum(prediccion_knn_cv == mhp_train$precio_bin) /nrow(mhp_train)
recall[i] = sum(prediccion_knn_cv == mhp_train$precio_bin & mhp_train$precio_bin == TRUE) / sum(mhp_train$precio_bin == TRUE)
precision[i] = sum(prediccion_knn_cv == mhp_train$precio_bin & prediccion_knn_cv == TRUE) / sum(prediccion_knn_cv == TRUE)
f1score[i] = 2*precision[i]*recall[i]/(precision[i]+recall[i])
}
resultados_knn = as.data.frame(cbind(accuracy,f1score,precision,recall))
resultados_knn = resultados_knn %>% mutate(index=as.factor(seq(1:long)))
max(resultados_knn$f1score)## [1] NaNwhich.max(resultados_knn$f1score)## integer(0)ggplot(data=resultados_knn,aes(x=index,y=accuracy)) +
geom_col(colour="cyan4",fill="cyan3")+
ggtitle("Accuracy")
ggplot(data=resultados_knn,aes(x=index,y=f1score)) +
geom_col(colour="orange4",fill="orange3") +
ggtitle("F1_score values")
Tratamos de encontrar el valor óptimo de k utilizando validación cruzada. Se calculan varias métricas de evaluación, como precisión, recall y F1-score, para cada valor de k y finalmente tratamos de graficar la precisión en función del valor de k, pero analizando el gráfico no observamos ninguna diferencia, así que cogemos 5, que es el tamaño por convención que se suele coger.
Pasamos a realizar predicciones en el conjunto de datos de entrenamiento y en el de test además de realizar la matriz de confusión.
# En train
prediccion_knn5_train =knn.cv(mhp_train[,c("exterior","rooms","m2","elevator", "garage")],
k=5, cl=mhp_train$precio_bin)
confusionMatrix(table(prediccion_knn5_train,mhp_train$precio_bin), positive= "1")## Confusion Matrix and Statistics
##
##
## prediccion_knn5_train 1 0
## 1 3233 531
## 0 630 3338
##
## Accuracy : 0.8498
## 95% CI : (0.8417, 0.8577)
## No Information Rate : 0.5004
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.6997
##
## Mcnemar's Test P-Value : 0.004026
##
## Sensitivity : 0.8369
## Specificity : 0.8628
## Pos Pred Value : 0.8589
## Neg Pred Value : 0.8412
## Prevalence : 0.4996
## Detection Rate : 0.4181
## Detection Prevalence : 0.4868
## Balanced Accuracy : 0.8498
##
## 'Positive' Class : 1
## #En test
prediccion_knn5_test=knn(mhp_train[,c("exterior","rooms","m2","elevator", "garage")], mhp_test[,c("exterior","rooms","m2","elevator", "garage")],
k=5, cl=mhp_train$precio_bin)
confusionMatrix(table(prediccion_knn5_test,mhp_test$precio_bin), positive= "1")## Confusion Matrix and Statistics
##
##
## prediccion_knn5_test 1 0
## 1 1577 250
## 0 356 1683
##
## Accuracy : 0.8432
## 95% CI : (0.8314, 0.8546)
## No Information Rate : 0.5
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.6865
##
## Mcnemar's Test P-Value : 1.996e-05
##
## Sensitivity : 0.8158
## Specificity : 0.8707
## Pos Pred Value : 0.8632
## Neg Pred Value : 0.8254
## Prevalence : 0.5000
## Detection Rate : 0.4079
## Detection Prevalence : 0.4726
## Balanced Accuracy : 0.8432
##
## 'Positive' Class : 1
## La matriz de confusión en train nos muestra la cantidad de predicciones correctas e incorrectas de cada clase (1 y 0) en el conjunto de entrenamiento.
Accuracy (exactitud): 0.8482. Es la proporción de predicciones correctas en el conjunto de entrenamiento. En otras palabras, el modelo KNN clasifica correctamente el 84.82% de las observaciones del conjunto de entrenamiento.
Kappa: 0.6963. El índice kappa mide la concordancia entre las predicciones y los valores reales. Un kappa de 1 indica una concordancia perfecta, mientras que un kappa de 0 indica una concordancia totalmente aleatoria. Nuestro kappa de 0.6963 muestra una buena concordancia entre las predicciones y los valores reales en el conjunto de entrenamiento.
Mcnemar’s Test P-Value: 0.003201. Esta prueba evalúa si hay una diferencia significativa entre las predicciones falsas positivas y las predicciones falsas negativas.
La matriz de confusión en test nos muestra valores muy parecidos en los indicadores anteriores lo que nos lleva a pensar que el modelo tiene un rendimiento consistente y generaliza bien lo datos ya que el modelo clasifica correctamente las observaciones y no parece sobreajustado.
Hemos realizado diferentes pruebas quitando variables y se observa que si quitamos la variable m2 nos dice que hay muchos empates. Sin embargo, quitando las otras variables el acierto en train y test varía.
Obteniendo más o menos un 84% de acierto creemos que nos puede servir para clasificar las nuevas casas que entrasen en el dataset, puesto que no hemos sido capaces de mejorar ese %.
10.3 Decision Trees
El modelo de Árbol de Decisión es un algoritmo de aprendizaje supervisado que funciona dividiendo iterativamente el conjunto de datos en subconjuntos basados en las características del mismo. Al final de este proceso, se genera un árbol con nodos y hojas que representan las divisiones y las decisiones.
Algunas características importantes de los árboles de decisión son la interpretabilidad, capacidad para manejar datos no lineales, el poco preprocesamiento de los datos, la sensibilidad a datos desequilibrados y el posible sobreajuste en el caso de hacer un arbol con demasiadas ramificaciones.
A continuación hacemos el modelo de arbol de decisión.
arbol <- rpart(precio_bin ~ garage + elevator + rooms + tipo_casa + precio_distrito, data = mhp_train, control = rpart.control(minsplit = 1))
fancyRpartPlot(arbol, sub = "")
Para evaluar el modelo utilizamos la función predict sobre el conjunto test usando el argumento type = “class” que indica que se espera una salida categórica creandose una tabla de contingencias “tab1” con las predicciones que coinciden o no respecto a las observaciones reales.
Calculamos y obtenemos el número 0.8396275 que representa la tasa de acierto (accuracy) del modelo. El Árbol de Decisión hizo predicciones correctas en casi el 84% de los casos en el conjunto de datos de prueba. Nos parece un porcentaje bastante aceptable.
tab1 = table(pred = predict(arbol, mhp_test, type = "class"),
obs = mhp_test$precio_bin)
ntest = nrow(mhp_test)
acierto1 = sum(diag(tab1))/ntest
acierto1## [1] 0.8396275Siguiendo con la validación del modelo podemos usar las funciones printcp y plotcp que se utilizan para analizar la complejidad y el rendimiento usando el “pruning” o poda de ramas del arbol que a menudo reduce su complejidad y permite una mejor generalización del mismo disminuyendo el subreajuste.
La función printcp es una función que muestra una tabla que contiene información sobre la complejidad y la tasa de error en diferentes niveles de tamaño del árbol proporcionando una visión general de cómo cambia el rendimiento a medida que aumenta la complejidad. El valor de cp ha de ser tal que la tasa de error de validación cruzada sea mínima.
La función plotcp crea un gráfico de la tasa de error. Este gráfico es útil para visualizar cómo se comporta el rendimiento del árbol a medida que aumenta su tamaño y ayuda a decidir qué tamaño de árbol es el más apropiado para el problema en cuestión.
printcp(arbol)##
## Classification tree:
## rpart(formula = precio_bin ~ garage + elevator + rooms + tipo_casa +
## precio_distrito, data = mhp_train, control = rpart.control(minsplit = 1))
##
## Variables actually used in tree construction:
## [1] elevator precio_distrito rooms
##
## Root node error: 3863/7732 = 0.49961
##
## n= 7732
##
## CP nsplit rel error xerror xstd
## 1 0.516438 0 1.00000 1.03055 0.0113762
## 2 0.044784 1 0.48356 0.48356 0.0097435
## 3 0.041677 3 0.39399 0.40098 0.0091108
## 4 0.016179 4 0.35232 0.35232 0.0086689
## 5 0.010000 6 0.31996 0.31996 0.0083418Obtenidos los resultados, pasamos a analizarlos:
- CP (Cost Complexity): Mide la complejidad asociado con cada nivel de poda en el árbol de decisión. Un CP más pequeño indica una complejidad menor.
- nsplit: Número de divisiones (splits) realizadas hasta el momento en el árbol.
- rel error: Error relativo en cada etapa de poda del árbol.
- xerror: Error de validación cruzada (cross-validated) estimado para el árbol en cada etapa de poda.
- xstd: Desviación estándar del error de validación cruzada (xerror).
En base a los resultados presentados, se puede observar que el árbol sin poda tiene un error relativo de 1.0 y un error de validación cruzada de 1.0189. A medida que se realizan más divisiones en el árbol el CP disminuye y también lo hace el error relativo y el error de validación cruzada.
Para seleccionar el árbol óptimo, generalmente se busca el CP más pequeño que tenga un error de validación cruzada dentro de una desviación estándar del mínimo error de validación cruzada.
arbol$cptable[which.min(arbol$cptable[, "xerror"]),
"CP"]## [1] 0.01El propósito es identificar el mejor árbol de decisión con respecto a su capacidad de generalización en datos no vistos. La función which.min() devuelve el índice del valor mínimo en la columna “xerror” de la tabla de complejidad cptable del árbol de decisión. Luego, con ese índice, se extrae el valor correspondiente de la columna “CP” en la tabla de complejidad.
plotcp(arbol)
Podamos el árbol como indica la función.
pruneTREE1 = prune(arbol, cp = arbol$cptable[which.min(arbol$cptable[,
"xerror"]), "CP"])
fancyRpartPlot(pruneTREE1, uniform = TRUE, main = "Pruned Classification Tree",
sub = "")
pruneTREE2 = prune(arbol, cp = arbol$cptable[5, "CP"])
fancyRpartPlot(pruneTREE2, uniform = TRUE, main = "Pruned Classification Tree",
sub = "")
tab2 = table(pred = predict(pruneTREE2, mhp_test, type = "class"),
obs = mhp_test$precio_bin)
acierto2 = sum(diag(tab2))/ntest
acierto2## [1] 0.8396275Tras la poda nos sale literalmente es mismo resultado.
Creamos la matriz de confusión.
prediccion_1 <- predict(arbol, newdata = mhp_train, type = "class")confusionMatrix(prediccion_1, mhp_train[["precio_bin"]])## Confusion Matrix and Statistics
##
## Reference
## Prediction 1 0
## 1 3246 619
## 0 617 3250
##
## Accuracy : 0.8401
## 95% CI : (0.8318, 0.8482)
## No Information Rate : 0.5004
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.6803
##
## Mcnemar's Test P-Value : 0.9773
##
## Sensitivity : 0.8403
## Specificity : 0.8400
## Pos Pred Value : 0.8398
## Neg Pred Value : 0.8404
## Prevalence : 0.4996
## Detection Rate : 0.4198
## Detection Prevalence : 0.4999
## Balanced Accuracy : 0.8401
##
## 'Positive' Class : 1
## Los resultados en la matriz de confusión son los siguientes:
- Verdaderos positivos (VP): 3246
- Verdaderos negativos (VN): 3250
- Falsos positivos (FP): 619
- Falsos negativos (FN): 617
También obtenemos:
- Accuracy (exactitud): 84.01% de predicciones correctas.
- Kappa: 68.03%
- Sensitivity (sensibilidad): 0.8403. También conocida como recall o tasa verdadera positiva, mide la proporción de casos positivos que se identificaron correctamente.
- Specificity (especificidad): 0.8400. Mide la proporción de casos negativos que se identificaron correctamente. Un valor más alto es mejor.
- Pos Pred Value (valor predictivo positivo): 0.8398. Mide la proporción de predicciones positivas que son realmente positivas. Un valor más alto es mejor.
- Neg Pred Value (valor predictivo negativo): 0.8404. Mide la proporción de predicciones negativas que son realmente negativas. Un valor más alto es mejor.
- Mcnemar’s Test P-Value: 0.9773. Esta prueba compara las diferencias entre los falsos positivos y los falsos negativos. Un valor alto de “p” (como en este caso) sugiere que no hay una diferencia significativa entre ellos.
En resumen, este modelo de árbol de decisión podado tiene una exactitud del 84.01% y un coeficiente kappa de 0.6803, lo que indica un buen rendimiento en la clasificación.
Sin embargo para este modelo se han tenido en cuenta bastantes variables, y aunque el rendimiento es bueno, se puede alcanzar un rendimiento casi igual utilizando sólo los metros cuadrados, como se puede ver a continuación
arbol <- rpart(precio_bin ~ m2, data = mhp_train, control = rpart.control(minsplit = 1))
fancyRpartPlot(arbol, sub = "")
tab1 = table(pred = predict(arbol, mhp_test, type = "class"),
obs = mhp_test$precio_bin)
ntest = nrow(mhp_test)
acierto1 = sum(diag(tab1))/ntest
acierto1## [1] 0.8344542Esta es una forma simple de confirmar la relevancia que tiene la superficie de una casa sobre su precio final, en contraste con otras variables como el tener ascensor o garaje. Además, esto nos deja un modelo mucho más simple y explicable.
10.4 Random Forest
Algoritmo de aprendizaje supervisado que se basa en la construcción de múltiples árboles de decisión durante el entrenamiento y la combinación de sus resultados para hacer una predicción. El objetivo principal de este enfoque es mejorar la precisión y la robustez del modelo en comparación con un único árbol de decisión.
El funcionamiento de un algoritmo Random Forest sería el siguiente:
Seleccionar muestras de arranque: El algoritmo selecciona aleatoriamente un subconjunto de muestras del conjunto de datos de entrenamiento, con reemplazo. Esto significa que una muestra puede ser seleccionada más de una vez.
Crear un árbol de decisión: Se crea un árbol de decisión utilizando el subconjunto de muestras seleccionado. Durante la construcción del árbol, en cada división del nodo, se selecciona un subconjunto aleatorio de características como candidatas para la división. La mejor división se elige en función de la medida de impureza como la entropía o el índice de Gini. Esto introduce aleatoriedad adicional en el proceso, lo que ayuda a reducir la correlación entre los árboles.
Repetir: El proceso se repite varias veces (número de árboles especificado) para construir un conjunto de árboles de decisión.
Combinar las predicciones: Para hacer una predicción en un nuevo caso, el algoritmo Random Forest ejecuta todos los árboles individualmente y luego combina sus resultados. En el caso de la clasificación, la clase con más votos (predicciones) de todos los árboles se elige como la predicción final.
Pasamos a construir el modelo comenzando con 10 árboles aleatórios.
# Escalamos la columna 4 en train y test
set.seed(19042023)
mhp_train_scaled <- mhp_train
mhp_train_scaled[, 4] <- scale(mhp_train[, 4])
mhp_test_scaled <- mhp_test
mhp_test_scaled[, 4] <- scale(mhp_test[, 4])
# Nos quedamos con todas las variables menos barrio y distrito porque hay muchas categorias
# Cogemos solo 10 árboles en el parámetro ntree
classifier <- randomForest(x = mhp_train_scaled[, c(3, 4, 5, 6, 7, 10,11)],
y = mhp_train_scaled$precio_bin,
ntree = 10)# Predicción de los resultados con el conjunto de testing
y_pred = predict(classifier, newdata = mhp_test_scaled[,c(3,4,5,6,7,10,11)])confusionMatrix(y_pred, mhp_test_scaled[["precio_bin"]])## Confusion Matrix and Statistics
##
## Reference
## Prediction 1 0
## 1 1676 127
## 0 257 1806
##
## Accuracy : 0.9007
## 95% CI : (0.8908, 0.9099)
## No Information Rate : 0.5
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.8013
##
## Mcnemar's Test P-Value : 4.61e-11
##
## Sensitivity : 0.8670
## Specificity : 0.9343
## Pos Pred Value : 0.9296
## Neg Pred Value : 0.8754
## Prevalence : 0.5000
## Detection Rate : 0.4335
## Detection Prevalence : 0.4664
## Balanced Accuracy : 0.9007
##
## 'Positive' Class : 1
## Ahora probamos con 100 árboles aleatorios para ver si cambia mucho el resultado.
set.seed(19042022)
classifier <- randomForest(x = mhp_train_scaled[, c(3, 4, 5, 6, 7, 10,11)],
y = mhp_train_scaled$precio_bin,
ntree = 100)y_pred = predict(classifier, newdata = mhp_test_scaled[,c(3,4,5,6,7,10,11)])confusionMatrix(y_pred, mhp_test_scaled[["precio_bin"]])## Confusion Matrix and Statistics
##
## Reference
## Prediction 1 0
## 1 1684 108
## 0 249 1825
##
## Accuracy : 0.9077
## 95% CI : (0.8981, 0.9166)
## No Information Rate : 0.5
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.8153
##
## Mcnemar's Test P-Value : 1.267e-13
##
## Sensitivity : 0.8712
## Specificity : 0.9441
## Pos Pred Value : 0.9397
## Neg Pred Value : 0.8799
## Prevalence : 0.5000
## Detection Rate : 0.4356
## Detection Prevalence : 0.4635
## Balanced Accuracy : 0.9077
##
## 'Positive' Class : 1
## La precisión del modelo es del 0,9077 y el intervalo de confianza del 95% es (0,8908, 0,9099), lo que indica que podemos estar bastante seguros de que la precisión real del modelo está dentro de este rango. El valor de kappa de 0,8013 indica que hay una buena concordancia entre las predicciones del modelo y las observaciones reales.
Además, la sensibilidad del modelo es del 0,8670 y la especificidad del 0,9343. La tasa de detección es del 0,4335 y la prevalencia es del 0,5000. En general, el modelo tiene un desempeño bastante bueno en la clasificación de las viviendas por encima y por debajo del umbral.
10.5 SVM
Support Vector Machine, algoritmo de aprendizaje supervisado que se utiliza tanto para problemas de clasificación como de regresión. En el caso de la clasificación, el objetivo principal de SVM es encontrar el hiperplano óptimo que separa los datos en diferentes clases. Un hiperplano es un límite de decisión que divide el espacio de características en dos partes, de modo que cada parte contenga puntos de datos de una clase específica. En el caso de la regresión, SVM busca encontrar una función que tenga el menor error de ajuste posible.
El algoritmo SVM funciona de la siguiente manera:
Encuentra el hiperplano óptimo: El algoritmo busca el hiperplano que tenga el margen máximo entre las dos clases. El margen es la distancia entre el hiperplano y los puntos de datos más cercanos a él, llamados “vectores de soporte”. Estos vectores de soporte son los puntos de datos que definen el límite de decisión y son fundamentales para el proceso de SVM.
Solución de un problema de optimización: La búsqueda del hiperplano óptimo implica la solución de un problema de optimización convexa. Este problema se resuelve mediante técnicas de optimización como el método de Lagrange o el algoritmo de optimización secuencial mínima (SMO).
Transformación a través del kernel: En casos en los que los datos no son linealmente separables en el espacio de características original, SVM puede utilizar una función de kernel para transformar los datos en un espacio de características de mayor dimensión donde se puedan separar linealmente. Los kernels comunes incluyen el kernel lineal, el kernel polinómico, el kernel de base radial (RBF) y el kernel sigmoide. La elección del kernel adecuado es crucial para el rendimiento del modelo SVM.
Clasificación: Para hacer una predicción en un nuevo punto de datos, el algoritmo SVM evalúa la posición del punto con respecto al hiperplano óptimo. Si el punto está por encima del hiperplano, se asigna a una clase, y si está por debajo, se asigna a la otra clase.
Para resolver de SVM necesito conocer el kernel y el parámetro de regularización que es una cuota superior sobre las variables duales, en definitiva, un coste. Pasamos a construir un primer el modelo usando un kernel radial RBF.
# Estandarizamos las variables numéricas "m2" y "rooms"
mhp_train$m2_estandarizado <- scale(mhp_train$m2)
mhp_train$rooms_estandarizado <- scale(mhp_train$rooms)
mhp_test$m2_estandarizado <- scale(mhp_test$m2)
mhp_test$rooms_estandarizado <- scale(mhp_test$rooms)
# Entrenamos el modelo SVM kernel con kernel radial
modelo_svm <- svm(precio_bin ~ rooms_estandarizado + m2_estandarizado + garage + elevator + tipo_casa + precio_distrito, data = mhp_train, kernel = "radial", cost = 10, scale = TRUE)
# Hacemos predicciones en el conjunto de prueba
svm_predicciones <- predict(modelo_svm, mhp_test[, c("rooms_estandarizado", "m2_estandarizado", "garage", "elevator", "tipo_casa", "precio_distrito")])
# Calculamos la matriz de confusión y la precisión del modelo
svm_matriz_confusion <- table(Prediccion = svm_predicciones, Real = mhp_test$precio_bin)
print(svm_matriz_confusion)## Real
## Prediccion 1 0
## 1 1692 116
## 0 241 1817precision <- sum(diag(svm_matriz_confusion)) / sum(svm_matriz_confusion)
print(paste0("Precisión del modelo SVM: ", round(precision * 100, 2), "%"))## [1] "Precisión del modelo SVM: 90.77%"Esto significa que:
- 1692 casas con precio alto fueron clasificadas correctamente como precio alto (verdaderos positivos).
- 1817 casas con precio bajo fueron clasificadas correctamente como precio bajo (verdaderos negativos).
- 241 casas con precio alto fueron clasificadas incorrectamente como precio bajo (falsos negativos).
- 116 casas con precio bajo fueron clasificadas incorrectamente como precio alto (falsos positivos).
La precisión del modelo es del 90.77%.
ELECCION DEL KERNEL
La elección del kernel en un modelo SVM depende de la naturaleza y la distribución de los datos. No hay una regla única para seleccionar el kernel, pero aquí hay algunas pautas generales para ayudarte a tomar una decisión:
Kernel lineal: Si tus datos son linealmente separables o la cantidad de atributos es alta en comparación con la cantidad de muestras, entonces un kernel lineal suele funcionar bien. Además, es computacionalmente más eficiente que otros kernels.
Kernel polinomial: Si tus datos tienen una estructura que puede ser capturada por un polinomio de un grado específico, el kernel polinomial podría ser una buena opción. Sin embargo, a medida que aumenta el grado del polinomio, el riesgo de sobreajuste aumenta y la complejidad computacional puede crecer.
Kernel radial (RBF): El kernel radial es una opción popular y versátil en muchos casos, especialmente cuando no se sabe de antemano si los datos son linealmente separables. Puede manejar datos no lineales y complejos. Sin embargo, encontrar los parámetros adecuados (costo y gamma) puede requerir una búsqueda en el espacio de parámetros y, en consecuencia, un mayor tiempo de cómputo.
Kernel sigmoide: El kernel sigmoide es similar a una función de activación de una red neuronal y puede ser una opción si deseas un enfoque similar al de una red neuronal pero con un modelo SVM.
Vamos a probar a representar los datos para ver si de alguna manera puede tomarse una mejor decisión sobre el Kernel a usar.
mhp_train %>%
ggplot(aes(x = m2, y = rooms, color = as.factor(precio_bin))) +
geom_point() +
facet_grid(elevator ~ garage, labeller = labeller(
elevator = c('0' = 'Sin Elevador', '1' = 'Con Elevador'),
garage = c('0' = 'Sin Garaje', '1' = 'Con Garaje')
)) +
labs(title = "Diagrama de dispersión de m2 vs rooms, según garage y elevator",
x = "Metros cuadrados",
y = "Número de habitaciones",
color = "Precio bin") +
scale_color_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"),
plot.title.position = "plot")
Probamos varios KERNELS
Creamos una función entrenar_evaluar_svm() que entrenará y evaluará un modelo SVM utilizando el tipo de kernel especificado en el argumento kernel_type. Luego con lapply() se aplicará esta función a cada tipo de kernel en el vector kernels. La función presentará un matriz de confusión y la precisión para cada tipo de kernel.
# Definir una función para entrenar y evaluar un modelo SVM con un tipo de kernel específico
entrenar_evaluar_svm <- function(kernel_type) {
# Entrenamos el modelo SVM con el kernel especificado
modelo_svm <- svm(precio_bin ~ rooms_estandarizado + m2_estandarizado + garage + elevator + tipo_casa + precio_distrito, data = mhp_train, kernel = kernel_type, cost = 10, scale = TRUE)
# Hacemos predicciones en el conjunto de prueba
svm_predicciones <- predict(modelo_svm, mhp_test[, c("rooms_estandarizado", "m2_estandarizado", "garage", "elevator", "tipo_casa", "precio_distrito")])
# Calculamos la matriz de confusión y la precisión del modelo
svm_matriz_confusion <- table(Prediccion = svm_predicciones, Real = mhp_test$precio_bin)
precision <- sum(diag(svm_matriz_confusion)) / sum(svm_matriz_confusion)
# Devolvemos una lista con la precisión y la matriz de confusión
return(list(precision = precision, matriz_confusion = svm_matriz_confusion))
}
# Probar diferentes tipos de kernel
kernels <- c("linear", "polynomial", "radial", "sigmoid")
resultados <- lapply(kernels, entrenar_evaluar_svm)
# Imprimir los resultados
for (i in 1:length(kernels)) {
cat(paste("Precisión del modelo SVM con kernel", kernels[i], ":", round(resultados[[i]]$precision * 100, 2), "%\n"))
print(resultados[[i]]$matriz_confusion)
}## Precisión del modelo SVM con kernel linear : 90.46 %
## Real
## Prediccion 1 0
## 1 1706 142
## 0 227 1791
## Precisión del modelo SVM con kernel polynomial : 90.43 %
## Real
## Prediccion 1 0
## 1 1688 125
## 0 245 1808
## Precisión del modelo SVM con kernel radial : 90.77 %
## Real
## Prediccion 1 0
## 1 1692 116
## 0 241 1817
## Precisión del modelo SVM con kernel sigmoid : 85 %
## Real
## Prediccion 1 0
## 1 1633 280
## 0 300 1653Se observa que el KERNEL RADIAL es el que mejor predice, ahora vamos a ajustar algunos parámetros del modelo utilizando la función tune para encontrar los mejores parámetros.
# Entrenar un modelo SVM con diferentes parámetros
modelo_svm_radial <- svm(precio_bin ~ rooms_estandarizado + m2_estandarizado + garage + elevator + tipo_casa + precio_distrito, data = mhp_train, kernel = "radial", cost = 10, scale = TRUE)
# Utilizar la función tune() para encontrar los mejores parámetros
tuned_parameters <- tune(svm,
precio_bin ~ rooms_estandarizado + m2_estandarizado + garage + elevator + tipo_casa + precio_distrito,
data = mhp_train,
kernel = "radial",
ranges = list(cost = c(0.1, 1, 10, 100, 1000),
gamma = c(0.1, 0.5, 1, 2, 5)))
# Entrenar un modelo SVM con los mejores parámetros encontrados
best_model <- tuned_parameters$best.model
# Hacer predicciones en el conjunto de prueba con el mejor modelo
best_predictions <- predict(best_model, mhp_test[, c("rooms_estandarizado", "m2_estandarizado", "garage", "elevator", "tipo_casa", "precio_distrito")])
# Calcular la matriz de confusión y la precisión del mejor modelo
best_matriz_confusion <- table(Prediccion = best_predictions, Real = mhp_test$precio_bin)
print(best_matriz_confusion)## Real
## Prediccion 1 0
## 1 1692 118
## 0 241 1815best_precision <- sum(diag(best_matriz_confusion)) / sum(best_matriz_confusion)
print(paste0("Precisión del mejor modelo SVM: ", round(best_precision * 100, 2), "%"))## [1] "Precisión del mejor modelo SVM: 90.71%"La precisión del modelo no mejora del 90.77% anterior y el coste computacional en la busqueda del mejor modelo es enorme. Parece que los parámetros del modelo elegidos durante el proceso de afinación no han sido los correctos. La afinación e hiperparámetros implica buscar la mejor combinación de hiperparámetros dentro del rango de búsqueda que se proporciona. Si el rango de búsqueda no contiene la combinación óptima de hiperparámetros, el proceso de afinación podría terminar eligiendo una combinación subóptima que resulte en un rendimiento inferior.
11 Evaluación y comparación de modelos
Tras la realización de todos los modelos, juntamos los mejores resultados en una sóla tabla para mejor comparabilidad: